text:
module TrainUtils where

import Control.Lens hiding (para)
import System.Random
import Control.Monad
import Prelude hiding (id, (.))
import Numeric.LinearAlgebra.Array
import Numeric.LinearAlgebra.Array.Util

import CategoricDefinitions
import Autodiff.GAD
import Autodiff.Dual
import Autodiff.Additive
import Autodiff.D
import Autodiff.Cont
import Ops
import Para
import TensorUtils
import OnesLike

-------------------------------------------
showNNInfo :: (ArrShow p, ArrShow a, ArrShow b)
    => Int -> a -> b -> LearnerType p a b -> IO ()
showNNInfo n a b nn = do
    putStrLn "-------------------------"
    putStrLn $ "Step " ++ arrShow n
    putStrLn $ "p\n" ++ arrShow (nn ^. p)
    putStrLn $ "a\n" ++ arrShow a
    putStrLn $ "a\n" ++ arrShow b

-- | Supervised learning training
-- Takes in a Learner, input-output pairs and a cost function
-- it partially applies the output to the cost function and composes the result inside learner
trainStepWithCost :: (OnesLike c, _)
    => LearnerType p a b -> (Int, IO (a, b), DType (b, b) c) -> IO (LearnerType p a b)
trainStepWithCost l (step, dataSampler, cost) = do
    (i, o) <- dataSampler
    when (step `mod` 100 == 0) $ showNNInfo step i o l
    let cost' = partiallyApply cost o
        (pGrad, _) = grad (cost' . (l ^. para . fn)) (l ^. p, i)
    return $ l & p .~ (l ^. optimizer) (l ^. p, pGrad)

instance (Random a, Random b) => Random (a, b) where
    random gen1 = let (x, gen2) = random gen1
                      (y, gen3) = random gen2
                  in ((x, y), gen3)

    randomR ((x1, y1), (x2, y2)) gen1 = let (x, gen2) = randomR (x1, x2) gen1
                                            (y, gen3) = randomR (y1, y2) gen2
                                        in ((x, y), gen3)


text:
{-# LANGUAGE FlexibleInstances #-}
{-# OPTIONS -Wall #-}
module NN.NeuralNetwork
  (
    NN (..)
  , NNVars
  , createNetworkW
  , Net (..)
  , module Util.Vars
  ) where

import AI.HNN.FF.Network
import Numeric.LinearAlgebra.HMatrix hiding (corr)
import System.Random (RandomGen, random)
import qualified Data.Vector as V
import Control.Monad.Reader
import RandomUtil.Random
import Util.Vars

import Matrix.Traversable
import Data.Traversable

type NNVars = Vars

class NN n where
  create :: NNVars -> IO n
  get :: n -> Vector Double -> Vector Double
  train :: n -> Sample Double -> n
  train net samp = trainV net [samp]
  trainV :: n -> Samples Double -> n
  randomize :: (RandomGen g) => g -> n -> (g,n)
  getWeights :: n -> V.Vector (Matrix Double)
  getVars :: n -> NNVars

data Net = Net
  {
    _network :: Network Double
  , _vars    :: NNVars
  }

instance NN Net where
  create vars = do
    n <- createNetworkW vars
    return $ Net n vars
  get net input = getNetwork (_vars net) (_network net) input
  trainV net io = let vars = _vars net in Net (trainNetwork vars (_network net) io) vars
  randomize g net = let vars = _vars net
                        (g',n) = randomizeNetwork vars g (_network net)
                    in (g', Net n vars)
  getWeights (Net (Network m) _) = m
  getVars (Net _ vars) = vars

createNetworkW :: NNVars -> IO (Network Double)
createNetworkW env = runReader createNetworkReader env

createNetworkReader :: Reader NNVars (IO (Network Double))
createNetworkReader = do
  numInput      <- asks (getVar numberInputsS)
  numHidden     <- asks (getVar numberHiddenS)
  numOutput     <- asks (getVar numberOutputsS)
  return $ createNetwork (round numInput) [(round numHidden)] (round numOutput)

getNetwork :: NNVars -> Network Double -> Vector Double -> Vector Double
getNetwork vars net input = runReader (getNetworkReader net input) vars

getNetworkReader :: Network Double -> Vector Double -> Reader NNVars (Vector Double)
getNetworkReader net input = do
  sigOrTanh <- asks (getVar sigmoidTanS)
  let (act, _) = getActivationFunctions sigOrTanh
  return $ output net act input

trainNetwork :: NNVars -> Network Double -> Samples Double -> Network Double
trainNetwork vars net samples = runReader (trainNetworkReader net samples) vars

trainNetworkReader :: Network Double -> Samples Double -> Reader NNVars (Network Double)
trainNetworkReader net samples = do
  timesToTrain  <- asks (getVar timesToTrainS)
  learningRate  <- asks (getVar learningRateS)
  sigOrTanh     <- asks (getVar sigmoidTanS)
  let (act, act') = getActivationFunctions sigOrTanh
  return $ trainNTimes (round timesToTrain) learningRate act act' net samples

type RandNet g = g -> Network Double -> (g, Network Double)

randomizeNetwork :: RandomGen g => NNVars -> RandNet g
randomizeNetwork vars g net = runReader (randomizeNetworkReader g net) vars

randomizeNetworkReader :: RandomGen g
  => g -> Network Double -> Reader NNVars (g, Network Double)
randomizeNetworkReader g net@(Network weights) = do
  mean <- asks (getVar meanS)
  std  <- asks (getVar stddevS)
  rate <- asks (getVar mutationRate)
  let mutF = probApply rate (addNormalNoise mean std)
      (g', weights') = mapAccumL (matrixTraverse mutF) g weights
  return $ (g',fromWeightMatrices weights')

type MutateElementF g a b = (g -> a -> (g, b))

-- With a probability of the provided double, apply the random function
-- Otherwise, apply nothing
probApply :: RandomGen g
  => Double -> MutateElementF g a a -> MutateElementF g a a
probApply p rF g a = let (r,g') = random g
                     in if r > p then randomID g' a else rF g' a

randomID :: MutateElementF g a a
randomID g a = (g,a)

getActivationFunctions :: (Floating a) => Double -> (ActivationFunction a, ActivationFunctionDerivative a)
getActivationFunctions inp = let asInt :: Integer
                                 asInt = round inp
                             in case asInt of
                                  1 -> (sigmoid, sigmoid')
                                  _ -> (tanh, tanh')

ext: hs

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
{-
 - *Main Data.Complex> h = concat [[1..5], replicate 11 0]
 - *Main Data.Complex> x = concat [[1..10], replicate 6 0]
 - *Main Data.Complex> convFft h x
 - [1.0,4.000000000000007,10.0,20.0,35.0,49.99999999999999,65.0,80.0,95.0,110.0,114.0,106.0,85.0,50.00000000000001,7.105427357601002e-15,7.105427357601002e-15]
 - *Main Data.Complex>
 -}
import Data.Complex

convOne :: (Num a) => [a] -> [a] -> Int -> a
convOne h x n    = sum (zipWith (*) kernel input)
    where kernel = drop (-begin) (reverse h)
          input  = drop begin x
          begin  = n - length h + 1

conv :: (Num a) => [a] -> [a] -> [a]
conv h x = [convOne h x n | n <- [0..(length x + length h - 2)]]

isPowerOfTwo :: Int -> Bool
isPowerOfTwo 1 = True
isPowerOfTwo x
    | mod x 2 == 0 = isPowerOfTwo (div x 2)
    | otherwise = False

fftRaw :: (RealFloat a) => [a] -> Int -> Int -> [Complex a]
fftRaw _ 0 _ = []
fftRaw [] _  _ = []
fftRaw (x0:_) 1 _ = [x0 :+ 0]
fftRaw x n s = zipWith (+) x1 x2 ++ (zipWith (-) x1 x2)
    where x1 = fftRaw x (div n 2) (2 * s)
          x2 = zipWith (*) [exp (0 :+ (-2 * pi * fromIntegral k / (fromIntegral n))) | k <- [0..((div n 2) - 1)]] (fftRaw (drop s x) (div n 2) (2 * s))

fft :: (RealFloat a) => [a] -> [Complex a]
fft x
    | isPowerOfTwo n = fftRaw x n 1
    | otherwise = error "FFT works only for powers of two"
    where n = length x

ifftRaw :: (RealFloat a) => [Complex a] -> Int -> Int -> [Complex a]
ifftRaw _ 0 _ = []
ifftRaw [] _  _ = []
ifftRaw (x0:_) 1 _ = [x0]
ifftRaw x n s = zipWith (+) x1 x2 ++ (zipWith (-) x1 x2)
    where x1 = ifftRaw x (div n 2) (2 * s)
          x2 = zipWith (*) [exp (0 :+ (2 * pi * fromIntegral k / (fromIntegral n))) | k <- [0..((div n 2) - 1)]] (ifftRaw (drop s x) (div n 2) (2 * s))

ifft :: (RealFloat a) => [Complex a] -> [a]
ifft x
    | isPowerOfTwo n = [realPart v / (fromIntegral n) | v <- ifftRaw x n 1]
    | otherwise = error "IFFT works only for powers of two"
    where n = length x

convFft :: (RealFloat a) => [a] -> [a] -> [a]
convFft h x
    | length x == length h = ifft (zipWith (*) (fft h) (fft x))
    | otherwise = error "Kernel and input data must have the same length"

max_stars_repo_name: sjappig/convhs

max_issues_repo_name: sjappig/convhs

max_forks_repo_name: sjappig/convhs


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}
import           Control.Monad
import           Control.Monad.Random
import           Data.List                    (foldl')

import qualified Data.ByteString              as B
import           Data.Semigroup               ((<>))
import           Data.Serialize

import           GHC.TypeLits

import qualified Numeric.LinearAlgebra.Static as SA

import           Options.Applicative

import           Grenade

-- The definition for our simple feed forward network.
-- The type level lists represents the layers and the shapes passed through the layers.
-- One can see that for this demonstration we are using relu, tanh and logit non-linear
-- units, which can be easily substitute for each other in and out.
--
-- With around 100000 examples, this should show two clear circles which have been learned by the network.
type FFNet = Network '[ FullyConnected 2 40, Tanh, FullyConnected 40 30, Relu, FullyConnected 30 20, Relu, FullyConnected 20 10, Relu, FullyConnected 10 1, Logit ]
                     '[ 'D1 2, 'D1 40, 'D1 40, 'D1 30, 'D1 30, 'D1 20, 'D1 20, 'D1 10, 'D1 10, 'D1 1, 'D1 1]

randomNet :: IO FFNet
randomNet = randomNetworkInitWith HeEtAl  -- you might want to try `Xavier` or `UniformInit` instead of `HeEtAl`

netTrain :: FFNet -> LearningParameters -> Int -> IO FFNet
netTrain net0 rate n = do
    inps <- replicateM n $ do
      s  <- getRandom
      return $ S1D $ SA.randomVector s SA.Uniform * 2 - 1
    let outs = flip map inps $ \(S1D v) ->
                 if v `inCircle` (fromRational 0.50, 0.50)  || v `inCircle` (fromRational (-0.50), 0.50)
                   then S1D $ fromRational 1
                   else S1D $ fromRational 0

    let trained = foldl' trainEach net0 (zip inps outs)
    return trained

  where trainEach !network (i,o) = train rate network i o

netLoad :: FilePath -> IO FFNet
netLoad modelPath = do
  modelData <- B.readFile modelPath
  either fail return $ runGet (get :: Get FFNet) modelData

renderClass :: IO ()
renderClass = do
  let testIns = [ [ (x,y)  | x <- [0..50] ]
                           | y <- [0..20] ]
  let outMat  = fmap (fmap (\(x,y) -> (render (x/25-1) (y/10-1)))) testIns
  putStrLn $ unlines outMat

  where
    render x y  | x == 0 && y == 0 = '+'
                | y == 0 = '-'
                | x == 0 = '|'
                | otherwise = let v = SA.vector [x,y] :: SA.R 2
                              in if v `inCircle` (fromRational 0.50, 0.50)  || v `inCircle` (fromRational (-0.50), 0.50)
                                 then '1'
                                 else ' '

netScore :: FFNet -> IO ()
netScore network = do
    let testIns = [ [ (x,y)  | x <- [0..50] ]
                             | y <- [0..20] ]
        outMat  = fmap (fmap (\(x,y) -> (render (x/25-1) (y/10-1) . normx) $ runNet network (S1D $ SA.vector [x / 25 - 1,y / 10 - 1]))) testIns
    putStrLn $ unlines outMat

  where
    render x y n'  | x == 0 && y == 0 = '+'
                   | y == 0 = '-'
                   | x == 0 = '|'
                   | n' <= 0.2  = ' '
                   | n' <= 0.4  = '.'
                   | n' <= 0.6  = '-'
                   | n' <= 0.8  = '='
                   | otherwise = '#'

normx :: S ('D1 1) -> Double
normx (S1D r) = SA.mean r

testValues :: FFNet -> IO ()
testValues network = do
  inps <- replicateM 1000 $ do
      s  <- getRandom
      return $ S1D $ SA.randomVector s SA.Uniform * 2 - 1
  let outs = flip map inps $ \(S1D v) ->
                 if v `inCircle` (fromRational 0.50, 0.50)  || v `inCircle` (fromRational (-0.50), 0.50)
                   then 1 :: Integer
                   else 0
  let ress = zip outs (map (round . normx . runNet network) inps)
      correct = length $ filter id $ map (uncurry (==)) ress
      incorrect = length $ filter id $ map (uncurry (/=)) ress
      falsePositives = length $ filter id $ map (uncurry (\shd nn -> shd == 0 && nn == 1)) ress
      falseNegatives = length $ filter id $ map (uncurry (\shd nn -> shd == 1 && nn == 0)) ress
  putStr $ show correct  ++ " | "
  putStr $ show incorrect ++ " | "
  putStr $ show falsePositives ++ " | "
  putStrLn $ show falseNegatives ++ " | "

inCircle :: KnownNat n => SA.R n -> (SA.R n, Double) -> Bool
v `inCircle` (o, r) = SA.norm_2 (v - o) <= r

data FeedForwardOpts = FeedForwardOpts Int LearningParameters

feedForward' :: Parser FeedForwardOpts
feedForward' =
  FeedForwardOpts <$> option auto (long "examples" <> short 'e' <> value 1000)
                  <*> (LearningParameters
                      <$> option auto (long "train_rate" <> short 'r' <> value 0.005)
                      <*> option auto (long "momentum" <> value 0.0)
                      <*> option auto (long "l2" <> value 0.0005)
                      )

main :: IO ()
main = do
  FeedForwardOpts examples rate <- execParser (info (feedForward' <**> helper) idm)

  putStrLn "| Nr | Correct | Incorrect | FalsePositives | FalseNegatives |"
  putStrLn "--------------------------------------------------------------"
  let nr = 100 :: Int
  mapM_ (\n -> do
    putStr $ "| " ++ show n  ++ " | "
    net0 <- randomNet
    net <- netTrain net0 rate examples
    -- netScore net
    testValues net) [1..nr]

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
module HNum.Special where

import           System.Random
import qualified Numeric.SpecFunctions         as SF
import           HNum.Vector
import           HNum.F

class FuncTools f => SpecialFunc f where
  erf :: Real a => f a -> f Double
  invErf :: Real a => f a -> f Double

instance SpecialFunc Vector where
  erf v = SF.erf . realToFrac <$> v
  invErf v = SF.invErf . realToFrac <$> v

max_stars_repo_name: Axect/HNumeric

max_issues_repo_name: Axect/HNumeric

max_forks_repo_name: Axect/HNumeric


text:
{-# LANGUAGE FlexibleContexts,
             UndecidableInstances,
             TypeFamilies #-}

module Network.Network
( Network(..) ) where

import Network.Layer
import System.Random
import Numeric.LinearAlgebra

-- | A network is
class Network a where
  predict :: Vector Double -> a -> Vector Double
  createNetwork :: (RandomGen g) => RandomTransform -> g -> [LayerDefinition] -> a

max_stars_repo_name: AkatsukiSirius/LambdaNet

max_issues_repo_name: world-admin/LambdaNet

max_forks_repo_name: world-admin/LambdaNet


text:
module Main where

import Control.Applicative
import Graphics.Rendering.Chart.Backend.Diagrams (toFile)
import qualified Graphics.Rendering.Chart.Easy as Plt
import Kalman
import Numeric.LinearAlgebra
import Sim

vanderPol = toFile Plt.def "VanderPol.svg" $ do
  let x0 = [1.4, 0]
  let q = diagl (replicate 2 0.0001)
  let r = diagl (replicate 2 0.5)
  let s0 = Estimate (vector [0, 5]) (diagl (replicate 2 5))

  let plant = Plant fA2 id fH2 id
  let sim = take 500 $ simulate plant (vector x0) 0
  let measurements = map snd sim
  let trajectory = map fst sim

  let sys = System fA2 fH2 q r
  let ekfs = scanl (ekf sys) s0 measurements
  let ukfs = scanl (ukf sys) s0 measurements

  let t1 = linspace 500 (0, 50 :: Double)
  let toPoints v = oneList (toList t1) $ map ((\[x] -> x) . toList . subVector 1 1) v

  Plt.plot (Plt.line "Plant" [toPoints trajectory])
  Plt.plot (Plt.line "UKF" [toPoints $ map sX ukfs])
  Plt.plot (Plt.line "EKF" [toPoints $ map sX ekfs])
  Plt.plot (Plt.line "Measurements" [toPoints measurements])

-- Plt.plot (Plt.line "LPF" [toPoints (movingAverage measurements)])

fA2 :: (Floating a) => [a] -> [a]
fA2 [x1, x2] = [x1 - x2 * dt, x2 + (- u * (1 - x1 ** 2) * x2 + x1) * dt]
  where
    u = 0.2
    dt = 0.1

fH2 :: (Floating a) => [a] -> [a]
fH2 [x1, x2] = [x1, x2]

oneList :: [a] -> [b] -> [(a, b)]
oneList [] _ = []
oneList (x : xs) (y : ys) = (x, y) : oneList xs ys

main = vanderPol


text:
{-# LANGUAGE BangPatterns #-}

-- |
-- Module    : Statistics.Matrix.Fast.Algorithms
-- Copyright : 2019 Magalame
-- License   : BSD3
--
-- Useful matrix functions.

module Statistics.Matrix.Fast.Algorithms
    (
      qr
    ) where

import Control.Applicative ((<$>), (<*>))
import Control.Monad.ST (ST, runST)
import Prelude hiding (replicate)
import Statistics.Matrix (Matrix (..),dimension, for)
import qualified Statistics.Matrix.Mutable as M
import qualified Data.Vector.Unboxed as U

-- | /O(r*c)/ Compute the QR decomposition of a matrix.
-- The result returned is the matrices (/q/,/r/).
qr :: Matrix -> (Matrix, Matrix)
qr mat = runST $ do
  let (m,n) = dimension mat

  r <- M.replicate n n 0
  a <- M.thaw mat
  for 0 n $ \j -> do
    cn <- M.immutably a $ \aa -> sqrt $ normCol j aa
    M.unsafeWrite r j j cn
    for 0 m $ \i -> M.unsafeModify a i j (/ cn)
    for (j+1) n $ \jj -> do
      p <- innerProduct a j jj
      M.unsafeWrite r j jj p
      for 0 m $ \i -> do
        aij <- M.unsafeRead a i j
        M.unsafeModify a i jj $ subtract (p * aij)
  (,) <$> M.unsafeFreeze a <*> M.unsafeFreeze r

normCol :: Int -> Matrix -> Double
normCol jthcol (Matrix r c v) = sub 0 0
  where sub !acc !ij | ij == r = acc
                     | otherwise = sub ( valCol*valCol + acc ) (ij+1)
                                   where
                                    valCol = U.unsafeIndex v (ij*c+jthcol)

innerProduct :: M.MMatrix s -> Int -> Int -> ST s Double
innerProduct mmat j k = M.immutably mmat $ \mat ->
  dotCol j mat k mat

dotCol :: Int -> Matrix -> Int -> Matrix -> Double
dotCol jthcol (Matrix r1 c1 v1) kthcol (Matrix _ c2 v2) = sub 0 0
  where sub !acc !ij | ij == r1 = acc
                     | otherwise = sub ( valColj*valColk + acc ) (ij+1)
                                   where
                                    valColk = U.unsafeIndex v2 (ij*c2+kthcol)
                                    valColj = U.unsafeIndex v1 (ij*c1+jthcol)

max_stars_repo_name: arvindd/dh-core

max_issues_repo_name: arvindd/dh-core

max_forks_repo_name: arvindd/dh-core


text:
{-# LANGUAGE CPP, TypeOperators, TypeFamilies #-}
{-# LANGUAGE DeriveFunctor, DeriveDataTypeable #-}
{-# OPTIONS_GHC -Wall #-}

-- {-# OPTIONS_GHC -fno-warn-unused-imports #-} -- TEMP
-- {-# OPTIONS_GHC -fno-warn-unused-binds   #-} -- TEMP

#define CustomComplex

#if !defined CustomComplex
{-# OPTIONS_GHC -fno-warn-orphans #-} -- for HasRep
#endif

----------------------------------------------------------------------
-- |
-- Module      :  ConCat.Complex
-- Copyright   :  (c) 2015 Conal Elliott and David Banas
-- License     :  BSD3
--
-- Maintainer  :  conal@conal.net
-- Stability   :  experimental
--
-- Simplified Complex type
----------------------------------------------------------------------

module ConCat.Complex (Complex(..),cis) where

#if defined CustomComplex
import Control.Applicative (liftA2)
import Data.Typeable
import Data.Data
import Test.QuickCheck (Arbitrary(..),CoArbitrary(..))
#else
import Data.Complex
#endif

#if defined CustomComplex
infixl 1 :+
data Complex a = a :+ a deriving (Functor,Eq,Show,Typeable,Data,Ord)

-- | The nonnegative magnitude of a complex number.
-- {-# SPECIALISE magnitude :: Complex Double -> Double #-}
magnitude :: (RealFloat a) => Complex a -> a

magnitude (x:+y) = sqrt (x*x + y*y)

-- From Data.Complex

-- magnitude (x:+y) =  scaleFloat k
--                      (sqrt (sqr (scaleFloat mk x) + sqr (scaleFloat mk y)))
--                     where k  = max (exponent x) (exponent y)
--                           mk = - k
--                           sqr z = z * z

-- | The phase of a complex number, in the range @(-'pi', 'pi']@.
-- If the magnitude is zero, then so is the phase.
-- {-# SPECIALISE phase :: Complex Double -> Double #-}
phase :: (RealFloat a) => Complex a -> a
phase (0 :+ 0)   = 0            -- SLPJ July 97 from John Peterson
phase (x:+y)     = atan2 y x

-- | Complex Identity Scaler
-- Returns a point on the unit circle, such that the angle between the
-- vector to that point and the positive x-axis is equal to the function
-- argument, given in radians.
cis :: (RealFloat a) => a -> Complex a
cis w = cos w :+ sin w

instance  (RealFloat a) => Num (Complex a)  where
    -- {-# SPECIALISE instance Num (Complex Float) #-}
    -- {-# SPECIALISE instance Num (Complex Double) #-}
    (x:+y) + (x':+y')   =  (x+x') :+ (y+y')
    (x:+y) - (x':+y')   =  (x-x') :+ (y-y')
    (x:+y) * (x':+y')   =  (x*x'-y*y') :+ (x*y'+y*x')
    negate (x:+y)       =  negate x :+ negate y
    abs z               =  magnitude z :+ 0
    signum (0:+0)       =  0
    signum z@(x:+y)     =  x/r :+ y/r  where r = magnitude z
    fromInteger n       =  fromInteger n :+ 0
    {-# INLINE (+)         #-}
    {-# INLINE (-)         #-}
    {-# INLINE (*)         #-}
    {-# INLINE negate      #-}
    {-# INLINE abs         #-}
    {-# INLINE signum      #-}
    {-# INLINE fromInteger #-}

instance  (RealFloat a) => Fractional (Complex a)  where
    -- {-# SPECIALISE instance Fractional (Complex Float) #-}
    -- {-# SPECIALISE instance Fractional (Complex Double) #-}
    recip z@(x:+y) = x/mag :+ (-y)/mag where mag = magnitude z
--     (x:+y) / (x':+y')   =  (x*x''+y*y'') / d :+ (y*x''-x*y'') / d
--                            where x'' = scaleFloat k x'
--                                  y'' = scaleFloat k y'
--                                  k   = - max (exponent x') (exponent y')
--                                  d   = x'*x'' + y'*y''

    fromRational a      =  fromRational a :+ 0

instance  (RealFloat a) => Floating (Complex a) where
    -- {-# SPECIALISE instance Floating (Complex Float) #-}
    -- {-# SPECIALISE instance Floating (Complex Double) #-}
    pi             =  pi :+ 0
    exp (x:+y)     =  expx * cos y :+ expx * sin y
                      where expx = exp x
    log z          =  log (magnitude z) :+ phase z

    sqrt (0:+0)    =  0
    sqrt z@(x:+y)  =  u :+ (if y < 0 then -v else v)
                      where (u,v) = if x < 0 then (v',u') else (u',v')
                            v'    = abs y / (u'*2)
                            u'    = sqrt ((magnitude z + abs x) / 2)

    sin (x:+y)     =  sin x * cosh y :+ cos x * sinh y
    cos (x:+y)     =  cos x * cosh y :+ (- sin x * sinh y)
    tan (x:+y)     =  (sinx*coshy:+cosx*sinhy)/(cosx*coshy:+(-sinx*sinhy))
                      where sinx  = sin x
                            cosx  = cos x
                            sinhy = sinh y
                            coshy = cosh y

    sinh (x:+y)    =  cos y * sinh x :+ sin  y * cosh x
    cosh (x:+y)    =  cos y * cosh x :+ sin y * sinh x
    tanh (x:+y)    =  (cosy*sinhx:+siny*coshx)/(cosy*coshx:+siny*sinhx)
                      where siny  = sin y
                            cosy  = cos y
                            sinhx = sinh x
                            coshx = cosh x

    asin z@(x:+y)  =  y':+(-x')
                      where  (x':+y') = log (((-y):+x) + sqrt (1 - z*z))
    acos z         =  y'':+(-x'')
                      where (x'':+y'') = log (z + ((-y'):+x'))
                            (x':+y')   = sqrt (1 - z*z)
    atan z@(x:+y)  =  y':+(-x')
                      where (x':+y') = log (((1-y):+x) / sqrt (1+z*z))

    asinh z        =  log (z + sqrt (1+z*z))
    acosh z        =  log (z + (z+1) * sqrt ((z-1)/(z+1)))
    atanh z        =  0.5 * log ((1.0+z) / (1.0-z))

instance Arbitrary a => Arbitrary (Complex a) where
  arbitrary = liftA2 (:+) arbitrary arbitrary
  shrink (x :+ y) = map (uncurry (:+)) (shrink (x,y))

instance CoArbitrary a => CoArbitrary (Complex a) where
  coarbitrary (x :+ y) = coarbitrary x . coarbitrary y

#endif

max_stars_repo_name: kenranunderscore/concat

max_issues_repo_name: con-kitty/concat

max_forks_repo_name: con-kitty/concat


text:
module Numeric.MixtureModel.Beta ( -- * General data types
                                   Sample
                                 , SampleIdx
                                 , Samples
                                 , ComponentIdx
                                 , Assignments
                                 , Weight
                                 -- * Beta parameters
                                 , BetaParam (..)
                                 , BetaParams
                                 , ComponentParams
                                 , paramFromMoments
                                 , paramToMoments
                                 , paramsFromAssignments
                                 , paramToMode
                                 -- * Beta distribution
                                 , Prob
                                 , betaProb
                                 -- * Gibbs sampling
                                 , estimateWeights
                                 , updateAssignments
                                 , updateAssignments'
                                 -- * Likelihood
                                 , likelihood
                                 -- * Classification
                                 , classify
                                 ) where

import           Data.Function (on)
import qualified Data.Vector as VB
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Unboxed as VU
import qualified Data.Vector.Unboxed.Mutable as MV
import           Control.Monad.ST

import           Numeric.Log hiding (sum)
import           Numeric.SpecFunctions (logBeta)
import           Statistics.Sample (meanVarianceUnb)

import           Data.Random
import           Data.Random.Distribution.Categorical

type Prob = Log Double
type Sample = Double
type SampleIdx = Int
type ComponentIdx = Int
type Weight = Double
data BetaParam = BetaParam { bpAlpha :: !Double -- ^ Alpha
                           , bpBeta  :: !Double -- ^ Beta
                           }
               deriving (Show, Eq, Ord)

-- k refers to number of components
-- N refers to number of samples
type Samples = VU.Vector Sample                         -- length == N
type Assignments = VU.Vector ComponentIdx               -- length == N
type BetaParams = VB.Vector BetaParam                   -- length == K
type ComponentParams = VB.Vector (Weight, BetaParam)    -- length == K

beta :: Double -> Double -> Log Double
beta a b = Exp $ logBeta a b

-- | `betaProb (a,b) x` is the probability of `x` under Beta
-- distribution defined by parameters `a` and `b`
betaProb :: BetaParam -> Sample -> Prob
betaProb (BetaParam a b) x =
    1/beta a b * realToFrac (x**(a-1)) * realToFrac ((1-x)**(b-1))

-- | Beta parameter from sample mean and variance
paramFromMoments :: Double            -- ^ mean
                 -> Double            -- ^ variance
                 -> Maybe BetaParam
paramFromMoments xbar v
  | c < 0      = Nothing
  | otherwise  = Just $ BetaParam (xbar * c) ((1 - xbar) * c)
  where c = xbar * (1 - xbar) / v - 1

-- | Mean and variance of the given beta parameter
paramToMoments :: BetaParam -> (Double, Double)
paramToMoments (BetaParam a b) =
  let mean = a / (a+b)
      var  = a*b / (a+b)^2 / (a+b+1)
  in (mean, var)

-- | The mode of the given beta parameter
paramToMode :: BetaParam -> Double
paramToMode (BetaParam a b) = (a - 1) / (a + b - 2)

-- | Beta parameter from samples
paramFromSamples :: VU.Vector Sample -> BetaParam
paramFromSamples v | V.null v = error "Numeric.MixtureModel.Beta: Can't estimate priors from no samples"
paramFromSamples v =
    case uncurry paramFromMoments $ meanVarianceUnb v of
      Just a  -> a
      Nothing -> error "Numeric.MixtureModel.Beta: Somehow we have a negative variance"

-- | Beta parameter for component given samples and their component assignments
paramFromAssignments :: Samples -> Assignments -> ComponentIdx -> BetaParam
paramFromAssignments samples assignments k =
  paramFromSamples $ V.map snd $ V.filter (\(k',_)->k==k') $ V.zip assignments samples

-- | Beta parameters for all components given samples and their component assignments
paramsFromAssignments :: Samples -> Int -> Assignments -> BetaParams
paramsFromAssignments samples ncomps assignments =
  V.fromList $ map (paramFromAssignments samples assignments) [0..ncomps-1]

-- | Draw a new assignment for a sample given beta parameters
drawAssignment :: ComponentParams -> Sample -> RVar ComponentIdx
drawAssignment params x =
  let probs = map (\(w,p)->realToFrac w * betaProb p x) $ V.toList params
  in case filter (isInfinite . ln . fst) $ zip probs [0..] of
        (x:_)     -> return $ snd x
        otherwise -> categorical
                     $ map (\(p,k)->(realToFrac $ p / sum probs :: Double, k))
                     $ zip probs [0..]

countIndices :: Int -> VU.Vector Int -> VU.Vector Int
countIndices n v = runST $ do
    accum <- V.thaw $ VU.replicate n 0
    V.forM_ v $ \k -> do n' <- MV.read accum k
                         MV.write accum k $! n'+1
    V.freeze accum

-- | Estimate the component weights of a given set of parameters
estimateWeights :: Assignments -> BetaParams -> ComponentParams
estimateWeights assignments params =
  let counts = countIndices (V.length params) assignments
      norm = realToFrac $ V.length assignments
      weights = V.map (\n->realToFrac n / norm) counts
  in V.zip (V.convert weights) params

-- | Sample assignments for samples under given parameters
updateAssignments' :: Samples -> ComponentParams -> RVar Assignments
updateAssignments' samples params =
  V.mapM (drawAssignment params) samples

-- | Gibbs update of sample assignments
updateAssignments :: Samples -> Int -> Assignments -> RVar Assignments
updateAssignments samples ncomps assignments =
  updateAssignments' samples
  $ estimateWeights assignments
  $ paramsFromAssignments samples ncomps assignments

-- | Likelihood of samples assignments under given model parameters
likelihood :: Samples -> ComponentParams -> Assignments -> Prob
likelihood samples params assignments =
    V.product ( V.map (\(k,x)->betaProb (snd $ params V.! k) x)
              $ V.zip assignments samples
              )
  * V.product (V.map (\k->realToFrac $ fst $ params V.! k) assignments)

-- | Maximum likelihood classification
classify :: ComponentParams -> Sample -> ComponentIdx
classify params x =
  fst
  $ V.maximumBy (compare `on` \(_,(w,p))->realToFrac w * betaProb p x)
  $ V.indexed params

max_stars_repo_name: bgamari/mixture-model

max_issues_repo_name: bgamari/mixture-model

max_forks_repo_name: bgamari/mixture-model


text:
-- Gaussian primes
-- http://www.codewars.com/kata/54d6abf84a35017d30000b26/

module Data.Complex.Gaussian.Prime where

import Control.Arrow ((&&&))
import Data.Complex.Gaussian (Gaussian (..), norm)

isGaussianPrime :: Gaussian -> Bool
isGaussianPrime g@(Gaussian a b) | a == 0 = f b
                                 | b == 0 = f a
                                 | otherwise = isPrime . norm $ g
                                 where f = uncurry (&&) . ((==3) . (`mod` 4) &&& isPrime) . abs
                                       isPrime n | n `elem` [2, 3, 5] = True
                                                 | n `mod` 2 == 0 || n `mod` 3 == 0 = False
                                                 | otherwise = notElem 0 . map (n `mod`) $ [5, 7 .. floor . sqrt . fromIntegral $ n]

max_stars_repo_name: gafiatulin/codewars

max_issues_repo_name: gafiatulin/codewars

max_forks_repo_name: gafiatulin/codewars


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators       #-}

module Grenade.Utils.Conversion
  ( -- Column major mode helper functions
    NrColumns
  , Row
  , Col
  , ColumnMajorModeIndex
  , columnMajorModeIndex
  , columnMajorModeRowCol

  , toLayerShape
  -- 1D conversions
  , toS1D
  , toS1DV
  , fromS1D
  , fromS1DV
  -- 2D converions
  , toS2D
  , toS2DV
  , fromS2D
  , fromS2DV
  -- further conversions
  , toRows
  , toCols
  , toColumnsS2D
  , toRowsS2D
  , fromRowMajorVectorToSD1
  , fromRowMajorVectorToSD1V
  , fromRowMajorVectorToSD2
  , fromRowMajorVectorToSD2V
  , fromColumnMajorVectorToSD2V
  ) where

import           Data.Maybe                   (fromJust, fromMaybe)
import           Data.Proxy
import           Data.Singletons.TypeLits
import qualified Data.Vector.Storable         as V
import           Foreign
import qualified Numeric.LinearAlgebra        as LA
import qualified Numeric.LinearAlgebra.Static as LAS
import           System.IO.Unsafe             (unsafePerformIO)
import           Unsafe.Coerce                (unsafeCoerce)

import           Grenade.Core.Shape
import           Grenade.Types
import           Grenade.Utils.Vector

-- Column major mode helper functions

type NrColumns = Int
type Row = Int
type Col = Int
type ColumnMajorModeIndex = Int

columnMajorModeIndex :: NrColumns -> Row -> Col -> ColumnMajorModeIndex
columnMajorModeIndex cols row col = row * cols + col

columnMajorModeRowCol :: NrColumns -> ColumnMajorModeIndex -> (Row, Col)
columnMajorModeRowCol cols idx = (idx `div` cols, idx `mod` cols)

-- import           Debug.Trace

-- test =
--   toCols $
--   (fromRowMajorVectorToSD2 (V.fromList [0..9]) :: S ('D2 2 5))
--   (fromRowMajorVectorToSD1 (V.fromList [0..9]) :: S ('D1 10))

-- testShapeV = fromRowMajorVectorToSD2V (V.fromList [0..9]) :: S ('D2 2 5)
-- testShape = fromRowMajorVectorToSD2 (V.fromList [0..9]) :: S ('D2 2 5)

-- | Convert the Shape to a list of rows.
toRows :: S x -> [V.Vector RealNum]
toRows x =
  case x of
    (S1DV v)  -> [v]
    (S1D v)   -> [unsafeCoerce v]
    (S2D m)   -> map LAS.extract . LAS.toRows $ m
    (S2DV {}) -> toRowsS2D x
    (S3D m)   -> map LAS.extract . LAS.toRows $ m
{-# INLINE toRows #-}

-- | Convert the Shape to a list of columns.
toCols :: S x -> [V.Vector RealNum]
toCols x =
  case x of
    S1DV v    -> map V.singleton . V.toList $ v
    S1D v     -> map V.singleton . V.toList $ (unsafeCoerce v :: V.Vector RealNum)
    (S2DV {}) -> toColumnsS2D x
    (S2D m)   -> map LAS.extract . LAS.toColumns $ m
    (S3D m)   -> map LAS.extract . LAS.toColumns $ m

-- | Converts the given vector to the correct layer shape.
toLayerShape :: S i -> S x -> S x
toLayerShape x y = case (x, y) of
  (S1D{}, S1DV{}) -> toS1D y
  (S1DV{}, S1D{}) -> fromS1D y
  (S2D{}, S2DV{}) -> toS2D y
  (S2DV{}, S2D{}) -> fromS2D y
  _               -> y
{-# INLINE toLayerShape #-}

-- 1D conversion

-- | Convert to S1D.
toS1D :: S ('D1 l) -> S ('D1 l)
toS1D (S1DV vec) = S1D $ (fromMaybe err . LAS.create $ vec)
  where
    err = error $ "wrong length of vector with " ++ show (V.length vec) ++ " in toS1D "
toS1D x@S1D{} = x
{-# INLINE toS1D #-}

-- | Convert from S1DV. This is the same as @toS1D@.
fromS1DV :: S ('D1 l) -> S ('D1 l)
fromS1DV = toS1D
{-# INLINE fromS1DV #-}

-- | Convert from S1D to S1DV.
fromS1D :: S ('D1 l) -> S ('D1 l)
fromS1D (S1D vec) = S1DV (LAS.extract vec)
fromS1D x@S1DV{}  = x
{-# INLINE fromS1D #-}

-- | Convert to S1DV (from S1D). Thsi is the same as @fromS1D@.
toS1DV :: S ('D1 l) -> S ('D1 l)
toS1DV = fromS1D
{-# INLINE toS1DV #-}

-- 2D conversions

-- | Convert from vector representation.
toS2D :: forall i j . S ('D2 i j) -> S ('D2 i j)
toS2D mat@S2DV{} = S2D $ LAS.matrix (concatMap V.toList $ toRows mat)
toS2D x@S2D{}    = x
{-# INLINE toS2D #-}

-- | Convert from S2DV. This is the same as @toS2D@.
fromS2DV :: S ('D2 i j) -> S ('D2 i j)
fromS2DV = toS2D
{-# INLINE fromS2DV #-}

-- | Convert from S2D to S2DV.
fromS2D :: S ('D2 i j) -> S ('D2 i j)
fromS2D (S2D mat) = S2DV . V.concat . map LAS.extract . LAS.toColumns $ mat
fromS2D x@S2DV{}  = x
{-# INLINE fromS2D #-}

-- | Convert to S2DV. This is the same as @fromS2D@.
toS2DV :: S ('D2 i j) -> S ('D2 i j)
toS2DV = fromS2D
{-# INLINE toS2DV #-}

-- test = LAS.create $ NLA.reshape 5 $ (V.fromList [0..9]) :: Maybe (LAS.L 2 5)
-- Just (matrix
--  [ 0.0, 1.0, 2.0, 3.0, 4.0
--  , 5.0, 6.0, 7.0, 8.0, 9.0 ] :: L 2 5)

-- | Efficiently extract the columns of the shape. Returns a list of j vectors.
toColumnsS2D :: forall i j . S ('D2 i j) -> [V.Vector RealNum]
toColumnsS2D (S2D mat) = map LAS.extract . LAS.toColumns $ mat
toColumnsS2D (S2DV vec) = map (\idx -> V.slice idx m vec) [0,m .. (V.length vec - m)]
  where
    m = fromIntegral $ natVal (Proxy :: Proxy i)
{-# INLINE toColumnsS2D #-}

-- | Efficiently extract the rows of the shape. Returns a list of i vectors.
toRowsS2D :: forall i j . S ('D2 i j) -> [V.Vector RealNum]
toRowsS2D (S2D mat) = map LAS.extract . LAS.toRows $ mat
toRowsS2D (S2DV vec) =
  unsafePerformIO $
  V.unsafeWith vec $ \from ->
    flip mapM [0 .. m - 1] $ \row -> do
      vec' <- createVector n
      V.unsafeWith vec' $ \to -> do
        let go (-1) = return ()
            go !col = do
              let idx = col * m + row
              x <- peekElemOff from idx
              pokeElemOff to col x
              go (col - 1)
        go (n - 1)
        return vec'
  where
    m = fromIntegral $ natVal (Proxy :: Proxy i)
    n = fromIntegral $ natVal (Proxy :: Proxy j)
{-# INLINE toRowsS2D #-}

-- | Convert from a row major vector to @SD1@.
fromRowMajorVectorToSD1 :: forall l . (KnownNat l) => V.Vector RealNum -> S ('D1 l)
fromRowMajorVectorToSD1 vec
  | V.length vec /= l = error $ "cannot create Vector R " ++ show l ++ " from vector with length " ++ show (V.length vec) ++ " in fromRowMajorVectorToSD1"
  | otherwise = S1D (unsafeCoerce vec)
  where l = fromIntegral $ natVal (Proxy :: Proxy l)
{-# INLINE fromRowMajorVectorToSD1 #-}

-- | Convert from a row major vector to @SD1V@.
fromRowMajorVectorToSD1V :: forall l . (KnownNat l) => V.Vector RealNum -> S ('D1 l)
fromRowMajorVectorToSD1V = S1DV
{-# INLINE fromRowMajorVectorToSD1V #-}

-- | Convert from a row major vector to @SD2@.
fromRowMajorVectorToSD2 :: forall i j . (KnownNat i, KnownNat j) => V.Vector RealNum -> S ('D2 i j)
fromRowMajorVectorToSD2 vec
  | V.length vec /= m * n = error $ "cannot create matrix L " ++ show (m,n) ++ " from vector length " ++ show (V.length vec) ++ " in fromRowMajorVectorToSD2"
  | otherwise = toS2D $ fromRowMajorVectorToSD2V vec
                -- in S2D $ unsafeCoerce v'
  where
    m = fromIntegral $ natVal (Proxy :: Proxy i)
    n = fromIntegral $ natVal (Proxy :: Proxy j)
{-# INLINE fromRowMajorVectorToSD2 #-}

-- toRowMajorVector :: (KnownNat i, KnownNat j) => S ('D2 i j) -> V.Vector RealNum
-- toRowMajorVector = V.concat . toColumnsS2D

-- | Convert from a row major vector to @SD2V@.
fromRowMajorVectorToSD2V :: forall i j . (KnownNat i, KnownNat j) => V.Vector RealNum -> S ('D2 i j)
fromRowMajorVectorToSD2V vec = unsafePerformIO $ do
  vec' <- createVector (V.length vec)
  V.unsafeWith vec $ \from ->
    V.unsafeWith vec' $ \to ->  do
    let go (-1) = return ()
        go !k = do
          -- let idx = nRow * m + nCol
              -- nCol = k `div` n
              -- nRow = k `mod` n
          let idx = columnMajorModeIndex n row col
              col = k - row * n
              row = k `mod` n
          x <- peekElemOff from k
          pokeElemOff to idx x
          go (k-1)
    go (V.length vec - 1)
  return $ S2DV vec'
  where
    m = fromIntegral $ natVal (Proxy :: Proxy i)
    n = fromIntegral $ natVal (Proxy :: Proxy j)
{-# INLINE fromRowMajorVectorToSD2V #-}

-- | Convert from a row major vector to @SD2V@.
fromColumnMajorVectorToSD2V :: forall i j . (KnownNat i, KnownNat j) => V.Vector RealNum -> S ('D2 i j)
fromColumnMajorVectorToSD2V = S2DV
{-# INLINE fromColumnMajorVectorToSD2V #-}

test =
  -- (\x -> fromColumnMajorVectorToSD2V x :: S ('D2 5 2)) $
  -- V.concat $
  -- toRowsS2D $
  -- ttoRowMajorVector $
  fromS2D $
  -- (\x -> fromRowMajorVectorToSD2V x :: S ('D2 5 2) ) $ V.concat $
  -- toRowsS2D $
  -- ((\(S2D x) -> S2DV $ V.concat $ map LAS.extract . LAS.toColumns $ x :: S ('D2 2 5)) )
  -- (fromRowMajorVectorToSD2V (V.fromList [0..9]) :: S ('D2 2 5))
  (fromRowMajorVectorToSD2 (V.fromList [0..9]) :: S ('D2 5 2))
  -- (fromRowMajorVectorToSD1 (V.fromList [0..9]) :: S ('D1 10))

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE FlexibleInstances, FlexibleContexts #-}
module Test.Gen where
import Test.Utils
import Test.QuickCheck
import Control.Monad
import qualified Numeric.LinearAlgebra as L
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector.Storable as V

squared_real_matrices :: Int -> Gen (L.Matrix Float)
squared_real_matrices k = do
  vs <- sequence (replicate (k*k) arbitrary) :: Gen [Float]
  return $ L.reshape k $ V.fromList vs

small_matrices :: Gen (L.Matrix Float)
small_matrices = do
  n <- choose (2,10)
  squared_real_matrices n

pair :: Gen a -> Gen b -> Gen (a,b)
pair = liftM2 (,)

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
module GalaxyStats
where

import qualified Data.Map as M
import Text.Printf

import Galaxy
import DataFunction
import Statistics
import Libaddutil.Named
import Utils

numStarSystemsInGalaxy :: Galaxy a -> Int
numStarSystemsInGalaxy = M.size . starsystems

numStarsInSystem :: StarSystem a -> Int
numStarsInSystem = M.size . stars

numStarsInGalaxy :: Galaxy a -> Int
numStarsInGalaxy g = sum $ map numStarsInSystem ((M.elems . starsystems) g)

starsPerStarSystemInGalaxy :: Galaxy a -> Float
starsPerStarSystemInGalaxy g = fromIntegral (numStarsInGalaxy g) / fromIntegral (numStarSystemsInGalaxy g)

bodiesPerStarsInGalaxy g = fromIntegral (numBodiesInGalaxy g) / fromIntegral (numStarsInGalaxy g)
planetsPerStarsInGalaxy g = fromIntegral (numPlanetsInGalaxy g) / fromIntegral (numStarsInGalaxy g)
bodiesPerStarsystemsInGalaxy g = fromIntegral (numBodiesInGalaxy g) / fromIntegral (numStarSystemsInGalaxy g)
planetsPerStarsystemsInGalaxy g = fromIntegral (numPlanetsInGalaxy g) / fromIntegral (numStarSystemsInGalaxy g)
moonsPerBodiesInGalaxy g = fromIntegral (numMoonsInGalaxy g) / fromIntegral (numBodiesInGalaxy g - numMoonsInGalaxy g)

numBodiesOrbitingStar :: Star a -> Int
numBodiesOrbitingStar s =
  let ps = (M.elems . planets) s
  in length ps + sum (map numSatellites ps)

numSatellites :: Planet a -> Int
numSatellites = M.size . satellites

numBodiesInStarSystem :: StarSystem a -> Int
numBodiesInStarSystem = sum . map numBodiesOrbitingStar . M.elems . stars

numBodiesInGalaxy :: Galaxy a -> Int
numBodiesInGalaxy = sum . map numBodiesInStarSystem . M.elems . starsystems

moonsInGalaxy :: Galaxy a -> [Planet a]
moonsInGalaxy = concatMap M.elems . map satellites . allBodies

numMoonsInGalaxy :: Galaxy a -> Int
numMoonsInGalaxy = length . moonsInGalaxy

numStarsByType t g = (filter (\s -> spectralType s == t)) (allStars g)

numBStarsInGalaxy = length . numStarsByType SpectralTypeB
numAStarsInGalaxy = length . numStarsByType SpectralTypeA
numFStarsInGalaxy = length . numStarsByType SpectralTypeF
numGStarsInGalaxy = length . numStarsByType SpectralTypeG
numKStarsInGalaxy = length . numStarsByType SpectralTypeK
numMStarsInGalaxy = length . numStarsByType SpectralTypeM

numBodiesByPlanetType a g = length $ (filter (\p -> planettype p == a)) (allBodies g)

numRockyBodiesByAtmosphere :: Atmosphere -> Galaxy a -> Int
numRockyBodiesByAtmosphere a g = length $ filter (\p -> planettype p == RockyPlanet a) (allBodies g)

isGasGiant p = case planettype p of
                 SmallGasGiant     -> True
                 MediumGasGiant    -> True
                 LargeGasGiant     -> True
                 VeryLargeGasGiant -> True
                 _                 -> False

numGasGiantsInGalaxy g = length $ (filter isGasGiant) (allBodies g)

numWaterWeatherSystemBodiesInGalaxy :: Galaxy a -> Int
numWaterWeatherSystemBodiesInGalaxy = numRockyBodiesByAtmosphere WaterWeatherSystem
numNoAtmosphereBodiesInGalaxy g = numBodiesByPlanetType NoAtmosphere g + numBodiesByPlanetType Planetoid g
numHabitableBodiesInGalaxy :: Galaxy a -> Int
numHabitableBodiesInGalaxy g = length $ filter sustainsLife (allBodies g)

showPerc :: Int -> Int -> String
showPerc a b = printf " (%.2f%%)" (100.0 * (fromIntegral a / fromIntegral b) :: Float)

showstars :: String -> Int -> Int -> String
showstars c st tot = c ++ show st ++ (showPerc st tot) ++ "\n"

planetMassesInGalaxy g = map planetMass (allBodies g)

minBodyMassInGalaxy = minimum . planetMassesInGalaxy

maxBodyMassInGalaxy = maximum . planetMassesInGalaxy

medBodyMassInGalaxy = median . planetMassesInGalaxy

avgBodyMassInGalaxy = average . planetMassesInGalaxy

planetTemperaturesInGalaxy g = map planetTemperature (allPlanets g)

minBodyTemperatureInGalaxy = kelvinToCelsius . minimum . planetTemperaturesInGalaxy

maxBodyTemperatureInGalaxy = kelvinToCelsius . maximum . planetTemperaturesInGalaxy

avgBodyTemperatureInGalaxy = kelvinToCelsius . averageInt . planetTemperaturesInGalaxy

medBodyTemperatureInGalaxy = kelvinToCelsius . median . planetTemperaturesInGalaxy

numPlanetsInGalaxy g = numBodiesInGalaxy g - numMoonsInGalaxy g

galaxyStats :: Galaxy a -> String
galaxyStats g =
  let starstotal = numStarsInGalaxy g
      bodiestotal = numBodiesInGalaxy g
      bstars = numBStarsInGalaxy g
      astars = numAStarsInGalaxy g
      fstars = numFStarsInGalaxy g
      gstars = numGStarsInGalaxy g
      kstars = numKStarsInGalaxy g
      mstars = numMStarsInGalaxy g
      gasbodies = numGasGiantsInGalaxy g
      waterbodies = numWaterWeatherSystemBodiesInGalaxy g
      noatmbodies = numNoAtmosphereBodiesInGalaxy g
      habitablebodies = numHabitableBodiesInGalaxy g
  in
     "Galaxy name: " ++ name g ++ "\n" ++
     "Number of star systems: " ++ (show . numStarSystemsInGalaxy) g ++ "\n" ++
     "Number of stars: " ++ (show . numStarsInGalaxy) g ++ "\n" ++
     "Number of stars / star system: " ++ (show2f . starsPerStarSystemInGalaxy) g ++ "\n" ++
     showstars "Number of stars of spectral type B: " bstars starstotal ++
     showstars "Number of stars of spectral type A: " astars starstotal ++
     showstars "Number of stars of spectral type F: " fstars starstotal ++
     showstars "Number of stars of spectral type G: " gstars starstotal ++
     showstars "Number of stars of spectral type K: " kstars starstotal ++
     showstars "Number of stars of spectral type M: " mstars starstotal ++
     "Number of bodies: " ++ (show . numBodiesInGalaxy) g ++ "\n" ++
     "Number of planets: " ++ (show . numPlanetsInGalaxy) g ++ "\n" ++
     "Number of planets / starsystem: " ++ (show2f . planetsPerStarsystemsInGalaxy) g ++ "\n" ++
     "Number of planets / star: " ++ (show2f . planetsPerStarsInGalaxy) g ++ "\n" ++
     "Number of moons: " ++ (show . numMoonsInGalaxy) g ++ "\n" ++
     "Number of moons / planet: " ++ (show2f . moonsPerBodiesInGalaxy) g ++ "\n" ++
     "Minimum body mass (Earths): " ++ (show . minBodyMassInGalaxy) g ++ "\n" ++
     "Median body mass (Earths): " ++ (show . medBodyMassInGalaxy) g ++ "\n" ++
     "Average body mass (Earths): " ++ (show . avgBodyMassInGalaxy) g ++ "\n" ++
     "Maximum body mass (Earths): " ++ (show . maxBodyMassInGalaxy) g ++ "\n" ++
     "Minimum body temperature (Celsius): " ++ (show . minBodyTemperatureInGalaxy) g ++ "\n" ++
     "Median body temperature (Celsius): " ++ (show . medBodyTemperatureInGalaxy) g ++ "\n" ++
     "Average body temperature (Celsius): " ++ (show . avgBodyTemperatureInGalaxy) g ++ "\n" ++
     "Maximum body temperature (Celsius): " ++ (show . maxBodyTemperatureInGalaxy) g ++ "\n" ++
     showstars "Number of gas giants: " gasbodies bodiestotal ++
     showstars "Number of oxygen bodies: " waterbodies bodiestotal ++
     showstars "Number of bodies without atmosphere: " noatmbodies bodiestotal ++
     showstars "Number of habitable bodies: " habitablebodies bodiestotal

max_stars_repo_name: anttisalonen/starrover

max_issues_repo_name: anttisalonen/starrover

max_forks_repo_name: anttisalonen/starrover


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE PartialTypeSignatures #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE OverloadedStrings #-}

module Model where

import           Control.Monad.Bayes.Class
import           Control.Monad.Bayes.Sampler
import qualified Data.Text.IO as TIO
import           Control.Monad
import           Control.Monad.Bayes.Inference.PMMH
import           Control.Monad.State
import Control.Monad.Bayes.Weighted
import           Numeric.Log
import           Text.Megaparsec

import           DataParser
import           Statistics
import           Utils

data Params = Params {
    rho :: Double, -- ^ Rate of detection
    beta :: Double, -- ^ Mean contact rate between susceptible and infected people
    gamma :: Double -- ^ Mean recovery rate
} deriving Show

data FixedParams = FixedParams {
    numPop :: Int,
    timeSlices :: Int
}

data LatentState = LatentState {
    sus :: Int, -- ^ Number of people susceptible to infection
    inf :: Int, -- ^ Number of people currently infected
    recov :: Int -- ^ Number of people recovered from infection
}

{-
observation model: Poisson rho * I
-}

-- | Model for how we observe the number of infected people
observationModel :: MonadSample m => Params -> LatentState -> m InfectionCount
observationModel (Params rho _ _) (LatentState _ inf _) =
    poisson (rho * fromIntegral inf)

{-
dN_SI <-
-}
-- | Transition the model a single time slice
transitionModelSingleStep
    :: MonadSample m => FixedParams -> Params -> LatentState -> m LatentState
transitionModelSingleStep (FixedParams numPop timeSlices) (Params rho beta gamma) (LatentState sus inf recov)
    = do
        let dt = 1 / fromIntegral timeSlices
        dN_SI <- binomial
            sus
            (1 - exp ((-beta * dt * (fromIntegral inf)) / (fromIntegral numPop))
            )
        dN_IR <- binomial
            inf
            (1 - exp (-gamma * dt)
            )
        let sus'   = sus - dN_SI
        let inf'   = inf + dN_SI - dN_IR
        let recov' = recov + dN_IR
        return (LatentState sus' inf' recov')

-- | Transition the model for a full step
transitionModel :: MonadSample m => FixedParams -> Params -> LatentState -> m LatentState
transitionModel fixedParams params =
    transitionModelSingleStep fixedParams params

-- | Simulate a single step, returning the new latent state and appending the observed infection count to the state.
simulateStep
    :: (MonadSample m, MonadState [InfectionCount] m)
    => FixedParams
    -> Params
    -> LatentState
    -> m LatentState
simulateStep fixedParams params latent = do
    latent'        <- transitionModel fixedParams params latent
    infectionCount <- observationModel params latent'
    modify (++ [infectionCount])
    return latent'

-- | Simulate nsteps steps of an epidemic with the specified initial state and parameters
simulateEpidemic
    :: MonadSample m => LatentState -> FixedParams -> Params -> Int -> m Epidemic
simulateEpidemic initialState fixedParams params nsteps =
    Epidemic
        <$> execStateT
                (repeatFunction nsteps (simulateStep fixedParams params) initialState)
                []

-- | Execute a single simulation of an epidemic
generateSingleEpidemic :: LatentState -> FixedParams -> Params -> Int -> IO Epidemic
generateSingleEpidemic initialState fixedParams params nsteps =
    sampleIO $ simulateEpidemic initialState fixedParams params nsteps

generateEpidemics :: LatentState -> FixedParams -> Params -> Int -> Int -> IO [Epidemic]
generateEpidemics initialState fixedParams params nsteps nepidemics = replicateM nepidemics (generateSingleEpidemic initialState fixedParams params nsteps)

fixedParams :: FixedParams
fixedParams = FixedParams 763 1

params :: Params
params = Params 0.9 2.0 0.6

initialState :: LatentState
initialState = LatentState 762 1 0

scoreEpidemicToDatum
    :: MonadInfer m
    => FixedParams
    -> Epidemic
    -> LatentState
    -> Params
    ->  m Params
scoreEpidemicToDatum fixedParams dat  x params = do
    let obs lambda y = score (poissonPdf lambda y)
        go [] x = return x
        go (y:ys) x = do
            x' <- transitionModel fixedParams params x
            obs ((fromIntegral $ inf x') * (rho params)) y
            return x'
    (go $! (unwrapEpidemic dat)) $! x
    return params

unwrapEpidemic :: Epidemic -> [Int]
unwrapEpidemic (Epidemic xs) = xs

{-
ddprior <- function(params) {
    dgamma(params[[1]], 0.3, 10, log = TRUE) +
    dgamma(params[[3]], 1, 8, log=TRUE) +
    dbeta(params[[2]], 2,7, log=TRUE)
}
-}
paramsPrior :: MonadSample m => m Params
paramsPrior = do
    pBeta <- Control.Monad.Bayes.Class.gamma 2 1
    pRho <- Control.Monad.Bayes.Class.beta 2.0 7.0
    pGamma <- Control.Monad.Bayes.Class.gamma 1.0  (1 / 8.0)
    return (Params pRho pBeta pGamma )

testInferenceEpidemic :: Int -> Int -> IO [[(Params, Numeric.Log.Log Double)]]
testInferenceEpidemic nsteps nparticles = do
    ys <- parseFromFile epidemicParser "data/datafile"
    case ys of (Left _) -> error "naughty"
               (Right dat) -> sampleIO $ do
                       pmmhRes <- prior $ pmmh nsteps 14 nparticles paramsPrior (scoreEpidemicToDatum fixedParams dat initialState)
                       return pmmhRes

parseFromFile p file = runParser p file <$> TIO.readFile file

extractParams :: (Params -> Double) -> [[(Params, Numeric.Log.Log Double)]] -> [Double]
extractParams project samples = (project . fst . head) <$> samples

--scoreEpidemicToData :: (MonadSample m, MonadState [InfectionCount] m => Epidemic ->  Params -> LatentState -> m Epidemic
--scoreEpidemicToData data params initialState = do
--    let obs lambda y = score (poissonPdf lambda y)
    -- simulate a new x from old x (x is a latent state)
    -- calculate new lambda
    -- score observation at time using lambda
    -- store x in the monad state
    -- repeat

generateSamples :: [Double]
generateSamples = sampleSTfixed $ replicateM 1000 (normal 3 1)

max_stars_repo_name: rossng/sir-monad

max_issues_repo_name: rossng/sir-monad

max_forks_repo_name: rossng/sir-monad


text:
{-# LANGUAGE RankNTypes   #-}
{-# LANGUAGE Strict       #-}
{-# LANGUAGE StrictData   #-}
module Pinwheel.BlockCudaMatrix where

import           Control.Concurrent.Async
import           Control.DeepSeq
import           Control.Monad            as M
import           Data.Array.Repa          as R
import           Data.Complex
import           Data.List                as L
import           Data.Vector.Storable     as VS
import           Data.Vector.Unboxed      as VU
import           Foreign.CUDA.BLAS        as BLAS
import           Foreign.CUDA.Driver      as CUDA
import           Text.Printf
import           Utils.BLAS
import           Utils.Parallel

data CuVec a
  = CuVecHost (VS.Vector a)
  | CuVecDevice (DevicePtr a)

data CuMat a =
  CuMat Int
        Int
        (CuVec a)

instance (Storable a) => Show (CuVec a) where
  show (CuVecHost vec) = printf "CuVecHost length %d" (VS.length vec)
  show (CuVecDevice _) = "CuVecDevice"

instance NFData (CuVec a) where
  rnf (CuVecHost vec) = seq vec ()
  rnf (CuVecDevice ptr) = seq ptr ()

instance (Storable a) => Show (CuMat a) where
  show (CuMat rows cols vec) =
    printf "CuMat rows = %d cols = %d %s" rows cols (show vec)

instance NFData (CuMat a) where
  rnf (CuMat rows cols vec) = rows `seq` cols `seq` vec `seq` ()

{-# INLINE getHostCuVec #-}
getHostCuVec :: (Storable a) => CuVec a -> VS.Vector a
getHostCuVec (CuVecHost vec) = vec
getHostCuVec (CuVecDevice _) = error "Error in getHostCuVec: no host vector."

{-# INLINE getHostCuMat #-}
getHostCuMat :: (Storable a) => CuMat a -> VS.Vector a
getHostCuMat (CuMat _ _ vec) = getHostCuVec vec

{-# INLINE getRowsCuMat #-}
getRowsCuMat :: CuMat a -> Int
getRowsCuMat (CuMat rows _ _) = rows

{-# INLINE getColsCuMat #-}
getColsCuMat :: CuMat a -> Int
getColsCuMat (CuMat _ cols _) = cols

{-# INLINE concatCuMat #-}
concatCuMat :: (Storable a) => [CuMat a] -> CuMat a
concatCuMat xs =
  let rows = L.foldl' (\s mat -> s + getRowsCuMat mat) 0 xs
      cols = getColsCuMat . L.head $ xs
  in if L.any (/= cols) . L.map getColsCuMat $ xs
       then error
              "concatCuMat: concat matries with different numbers of columns."
       else CuMat rows cols . CuVecHost . VS.concat . L.map getHostCuMat $ xs

{-# INLINE unsafeWithCuMat #-}
unsafeWithCuMat :: (Storable e) => CuMat e -> (DevicePtr e -> IO a) -> IO a
unsafeWithCuMat (CuMat _ _ (CuVecHost vec)) f = unsafeWithGPU vec f
unsafeWithCuMat (CuMat _ _ (CuVecDevice vec)) f = f vec

{-# INLINE freeCuMat #-}
freeCuMat :: CuMat e -> IO ()
freeCuMat (CuMat _ _ (CuVecHost _)) = return ()
freeCuMat (CuMat _ _ (CuVecDevice ptr)) = CUDA.free ptr

{-# INLINE subMatMul #-}
subMatMul ::
     (Storable e, Floating e, CUBLAS e)
  => Handle
  -> Int
  -> [CuMat e]
  -> CuMat e
  -> IO (CuMat e)
subMatMul handle rows matA@((CuMat rowsA colsA _):_) b@(CuMat _ colsB _) =
  fmap (CuMat rows colsB . CuVecHost . VS.concat) . unsafeWithCuMat b $ \bPtr ->
    M.mapM
      (\harmonic ->
         unsafeWithCuMat harmonic $ \harmonicPtr ->
           gemmCuBLAS handle rowsA colsB colsA harmonicPtr bPtr)
      matA

-- Given row-major Mat1: rows1xcols1 and Mat2: rows2xcols2, this function outputs Mat1 X Mat2
matMul ::
     (Storable e, Floating e, CUBLAS e, Unbox e)
  => Bool
  -> Int
  -> Int
  -> [CuMat e]
  -> [CuMat e]
  -> IO (CuMat e)
matMul doesTranspose deviceID rows matA' matB = do
  dev <- device deviceID
  ctx <- CUDA.create dev []
  handle <- BLAS.create
  matA <-
    if L.length matA' == 1
      then do
        let (CuMat rowsA colsA (CuVecHost vec)) = L.head matA'
            len = rowsA * colsA
        devPtr <- CUDA.mallocArray len
        unsafeWith vec $ \ptr -> CUDA.pokeArray len ptr devPtr
        return [CuMat rowsA colsA (CuVecDevice devPtr)]
      else return matA'
  coefs <- M.mapM (subMatMul handle rows matA) matB
  let transposedCoefs = concatCuMat . parMap rdeepseq transposeCuMat $ coefs
      output =
        if doesTranspose
          then transposedCoefs
          else transposeCuMat transposedCoefs
  when (L.length matA' == 1) (freeCuMat . L.head $ matA)
  BLAS.destroy handle
  CUDA.destroy ctx
  return output

-- matA is divided into [rowsA1 x colsA .. rowsAN x colsA]
-- matB is divided into [colsA x colsB1 .. colsA x colsBM]
-- each colsA x colsBm is further divided into [colsA x colsBmk]
-- M is the number of GPUs
{-# INLINE blockMatrixMultiply #-}
blockMatrixMultiply ::
     (Storable e, Floating e, CUBLAS e, Unbox e)
  => Bool
  -> [Int]
  -> [CuMat e]
  -> [[CuMat e]]
  -> IO (CuMat e)
blockMatrixMultiply doesTranspose deviceIDs matAs matBss = do
  unless
    (L.length deviceIDs == L.length matBss)
    (error $
     printf
       "Error in blockMatrixMultiply: matrix B is not divided according to the number of GPUs.\n%d GPUs vs %d blocks"
       (L.length deviceIDs)
       (L.length matBss))
  let rows = L.sum . L.map getRowsCuMat $ matAs
  print matAs
  print matBss
  initialise []
  output <-
    fmap concatCuMat .
    mapConcurrently
      (\(deviceID, matBs) -> matMul True deviceID rows matAs matBs) .
    L.zip deviceIDs $
    matBss
  return $!
    if doesTranspose
      then output
      else transposeCuMat output

-- Utilities
{-# INLINE transposeCuMat #-}
transposeCuMat :: (Storable e, Unbox e) => CuMat e -> CuMat e
transposeCuMat (CuMat rows cols (CuVecHost vec)) =
  CuMat cols rows .
  CuVecHost .
  VS.convert .
  toUnboxed .
  computeS .
  R.backpermute (Z :. cols :. rows) (\(Z :. r :. c) -> (Z :. c :. r)) .
  fromUnboxed (Z :. rows :. cols) . VS.convert $
  vec
transposeCuMat _ = error "transposeCuMat: Cannot transpose matrices on a device."

{-# INLINE createCuMat #-}
createCuMat :: (Storable e) => [VS.Vector e] -> CuMat e
createCuMat xs =
  CuMat (L.length xs) (VS.length . L.head $ xs) . CuVecHost . VS.concat $ xs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE BangPatterns #-}
module Polestar.Type where
import Data.Complex
import Data.Monoid
import Data.Ord

-- variable name (for debugging)
newtype Id = Id String deriving (Show)

-- ordering: s <: t implies s < t
data PrimType
  -- numeric types
  = PTyZero
  | PTyNat
  | PTyInt
  | PTyNNReal -- non-negative real number
  | PTyReal
  | PTyImaginary
  | PTyComplex
  -- boolean
  | PTyTrue
  | PTyFalse
  | PTyBool
  -- other
  | PTyUnit
  deriving (Eq,Ord,Show,Enum,Bounded)

data Type = TyPrim !PrimType
          | TyArr Type Type
          | TyAll Id !(Maybe Type) Type
          | TyRef !Int
          | TyTuple [Type] -- must have 2 or more elements
          | TyInter [Type] -- must have 2 or more elements
          deriving (Eq,Show)

pattern TyZero = TyPrim PTyZero
pattern TyNat = TyPrim PTyNat
pattern TyInt = TyPrim PTyInt
pattern TyNNReal = TyPrim PTyNNReal
pattern TyReal = TyPrim PTyReal
pattern TyImaginary = TyPrim PTyImaginary
pattern TyComplex = TyPrim PTyComplex
pattern TyTrue = TyPrim PTyTrue
pattern TyFalse = TyPrim PTyFalse
pattern TyBool = TyPrim PTyBool
pattern TyUnit = TyPrim PTyUnit
pattern TyPair a b = TyTuple [a,b]

-- compound canonical type
-- no pair of the components have subtyping relation
type CanonicalType = [ICanonicalType]

-- individual canonical type
data ICanonicalType = CTyPrim !PrimType
                    | CTyArr CanonicalType ICanonicalType
                    | CTyAll Id CanonicalType ICanonicalType
                    | CTyRef !Int
                    | CTyTuple [ICanonicalType]
                    deriving (Eq,Show)

pattern CTyTrue = CTyPrim PTyTrue
pattern CTyFalse = CTyPrim PTyFalse
pattern CTyBool = CTyPrim PTyBool

data Builtin
  -- unary
  = BNegate
  | BLogicalNot
  | BNatToInt
  | BNatToNNReal
  | BIntToNat -- max(n,0)
  | BIntToReal
  | BNNRealToReal
  | BRealToComplex
  | BMkImaginary
  | BImaginaryToComplex
  | BRealPart
  | BImagPart
  | BAbs
  | BSqrt
  | BExp
  | BExpm1
  | BLog
  | BLog1p
  | BSin
  | BCos
  | BTan
  | BSinh
  | BCosh
  | BTanh
  | BAsin
  | BAcos
  | BAtan
  | BAsinh
  | BAcosh
  | BAtanh
  -- TODO: factorial
  -- TODO: exp, log, expm1, log1p
  -- TODO: trigonometric functions and hyperbolic functions
  -- binary
  | BAdd
  | BSub
  | BMul
  | BDiv
  | BPow
  | BTSubNat -- truncated subtraction
  | BLt
  | BLe
  | BEqual
  | BMax
  | BMin
  | BIntDiv
  | BIntMod
  | BGcd
  | BLcm
  | BLogicalAnd
  | BLogicalOr
  -- TODO: binomial coefficients
  -- other primitives
  | BIterate
  -- | BUnsafeGlue
  deriving (Eq,Show,Enum,Bounded)

data PrimValue = PVZero
               | PVInt !Integer
               | PVReal !Double
               | PVImaginary !Double
               | PVComplex !(Complex Double)
               | PVBool !Bool
               | PVUnit
               | PVBuiltin !Builtin
               deriving (Eq,Show)

data Term = TmPrim !PrimValue             -- primitive value
          | TmAbs Id Type Term            -- lambda abstraction
          | TmTyAbs Id !(Maybe Type) Term -- bounded type abstraction
          | TmRef !Int                    -- variable (de Bruijn index)
          | TmApp Term Term               -- function application
          | TmTyApp Term Type             -- type application
          | TmLet Id Term Term            -- let-in
          | TmAlt Id [Type] Term          -- type alternation
          | TmIf Term Term Term           -- if-then-else
          | TmTuple [Term]                -- tuple
          | TmProj Term !Int              -- projection
          | TmCoerce Term Type            -- coercion
          | TmCoherentTuple [Term]        -- coherent tuple
          deriving (Eq,Show)

data Binding = VarBind Id CanonicalType
             | TyVarBind Id CanonicalType
             | AnonymousBind
             deriving (Eq,Show)

isUnary :: Builtin -> Bool
isUnary f = case f of
  BNegate -> True
  BLogicalNot -> True
  BNatToInt -> True
  BNatToNNReal -> True
  BIntToNat -> True
  BIntToReal -> True
  BNNRealToReal -> True
  BRealToComplex -> True
  BMkImaginary -> True
  BImaginaryToComplex -> True
  BRealPart -> True
  BImagPart -> True
  BAbs -> True
  BSqrt -> True
  BExp -> True
  BExpm1 -> True
  BLog -> True
  BLog1p -> True
  BSin -> True
  BCos -> True
  BTan -> True
  BSinh -> True
  BCosh -> True
  BTanh -> True
  BAsin -> True
  BAcos -> True
  BAtan -> True
  BAsinh -> True
  BAcosh -> True
  BAtanh -> True
  _ -> False

isBinary :: Builtin -> Bool
isBinary f = case f of
  BAdd -> True
  BSub -> True
  BMul -> True
  BDiv -> True
  BPow -> True
  BTSubNat -> True
  BLt -> True
  BLe -> True
  BEqual -> True
  BMax -> True
  BMin -> True
  BIntDiv -> True
  BIntMod -> True
  BGcd -> True
  BLcm -> True
  BLogicalAnd -> True
  BLogicalOr -> True
  _ -> False

isValue :: Term -> Bool
isValue t = case t of
  TmPrim _ -> True
  TmAbs _ _ _ -> True
  TmTyAbs _ _ _ -> True
  TmApp (TmPrim (PVBuiltin f)) x | isBinary f -> isValue x -- partial application
  TmTyApp (TmPrim (PVBuiltin BIterate)) ty -> True
  TmApp (TmTyApp (TmPrim (PVBuiltin BIterate)) ty) x -> isValue x
  TmApp (TmApp (TmTyApp (TmPrim (PVBuiltin BIterate)) ty) x) y -> isValue x && isValue y
  TmTuple xs -> all isValue xs
  TmCoherentTuple xs -> all isValue xs
  _ -> False

getCTypeFromContext :: [Binding] -> Int -> CanonicalType
getCTypeFromContext ctx i
  | 0 <= i && i < length ctx = case ctx !! i of
      VarBind _ ty -> ty
      b -> error ("TmRef: expected a variable binding, found " ++ show b)
  | otherwise = error "TmRef: index out of bounds"

getTypeFromContext :: [Binding] -> Int -> Type
getTypeFromContext ctx i = canonicalToOrdinary $ getCTypeFromContext ctx i

getCBoundFromContext :: [Binding] -> Int -> CanonicalType
getCBoundFromContext ctx i
  | i < length ctx = case ctx !! i of
                       TyVarBind _ b -> b
                       b -> error ("TyRef: expected a type variable binding, found " ++ show b)
  | otherwise = error "TyRef: index out of bounds"

getBoundFromContext :: [Binding] -> Int -> Maybe Type
getBoundFromContext ctx i = canonicalToOrdinaryM $ getCBoundFromContext ctx i

typeShift :: Int -> Int -> Type -> Type
-- typeShift 0 _ ty = ty
typeShift !delta = go
  where
    go !i ty = case ty of
      TyPrim _ -> ty
      TyArr u v -> TyArr (go i u) (go (i + 1) v)
      TyAll name b t -> TyAll name (typeShift delta i <$> b) (typeShift delta (i + 1) t)
      TyRef j | j >= i, j + delta >= 0 -> TyRef (j + delta)
              | j >= i, j + delta < 0 -> error "typeShift: negative index"
              | otherwise -> ty
      TyTuple tys -> TyTuple $ map (go i) tys
      TyInter tys -> TyInter $ map (go i) tys

typeShiftC :: Int -> Int -> CanonicalType -> CanonicalType
-- typeShiftC 0 _ = id
typeShiftC !delta !i = map (typeShiftI delta i)

typeShiftI :: Int -> Int -> ICanonicalType -> ICanonicalType
typeShiftI !delta = go
  where
    go :: Int -> ICanonicalType -> ICanonicalType
    go !i ty = case ty of
      CTyPrim _ -> ty
      CTyArr u v -> CTyArr (map (go i) u) (go (i + 1) v)
      CTyAll name b t -> CTyAll name (map (go i) b) (go (i + 1) t)
      CTyRef j | j >= i, j + delta >= 0 -> CTyRef (j + delta)
               | j >= i, j + delta < 0 -> error "typeShift: negative index"
               | otherwise -> ty
      CTyTuple tys -> CTyTuple $ map (go i) tys

typeSubstD :: Int -> Type -> Int -> Type -> Type
typeSubstD !depth s !i ty = case ty of
  TyPrim _ -> ty
  TyArr u v -> TyArr (typeSubstD depth s i u) (typeSubstD (depth + 1) s (i + 1) v)
  TyAll name b t -> TyAll name (typeSubstD depth s i <$> b) (typeSubstD (depth + 1) s (i + 1) t)
  TyRef j | j == i -> typeShift depth 0 s
          | j > i -> TyRef (j - 1)
          | otherwise -> ty
  TyTuple tys -> TyTuple $ map (typeSubstD depth s i) tys
  TyInter tys -> TyInter $ map (typeSubstD depth s i) tys

typeSubst = typeSubstD 0

typeSubstCD :: Int -> CanonicalType -> Int -> ICanonicalType -> CanonicalType
typeSubstCD !depth s !i ty = case ty of
  CTyPrim _ -> return ty
  CTyArr u v -> CTyArr (u >>= typeSubstCD depth s i) <$> typeSubstCD (depth + 1) s (i + 1) v
  CTyAll name b t -> CTyAll name (b >>= typeSubstCD depth s i) <$> typeSubstCD (depth + 1) s (i + 1) t
  CTyRef j | j == i -> typeShiftC depth 0 s
           | j > i -> return $ CTyRef (j - 1)
           | otherwise -> return ty
  CTyTuple tys -> CTyTuple <$> mapM (typeSubstCD depth s i) tys

typeSubstC = typeSubstCD 0

canonicalToOrdinary :: CanonicalType -> Type
canonicalToOrdinary [] = error "canonicalToOrdinary: Top type"
canonicalToOrdinary [t] = iCanonicalToOrdinary t
canonicalToOrdinary tys = TyInter $ map iCanonicalToOrdinary tys

canonicalToOrdinaryM :: CanonicalType -> Maybe Type
canonicalToOrdinaryM [] = Nothing
canonicalToOrdinaryM [t] = Just $ iCanonicalToOrdinary t
canonicalToOrdinaryM tys = Just $ TyInter $ map iCanonicalToOrdinary tys

iCanonicalToOrdinary :: ICanonicalType -> Type
iCanonicalToOrdinary (CTyPrim p) = TyPrim p
iCanonicalToOrdinary (CTyArr u v) = TyArr (canonicalToOrdinary u) (iCanonicalToOrdinary v)
iCanonicalToOrdinary (CTyAll name b t) = TyAll name (canonicalToOrdinaryM b) (iCanonicalToOrdinary t)
iCanonicalToOrdinary (CTyRef i) = TyRef i
iCanonicalToOrdinary (CTyTuple tys) = TyTuple $ map iCanonicalToOrdinary tys

instance Eq Id where
  _ == _ = True

{-
instance Ord ICanonicalType where
  compare (CTyPrim p) (CTyPrim p') = compare p p'
  compare (CTyArr s t) (CTyArr s' t') = compare (Down s) (Down s') <> compare t t'
  compare (CTyAll _ s t) (CTyAll _ s' t') = compare (Down s) (Down s') <> compare t t'
  compare (CTyRef i) (CTyRef i') = compare i i'
  compare (CTyTuple tys) (CTyTuple tys') = compare tys tys'
  compare (CTyPrim _) _ = LT
  compare _ (CTyPrim _) = GT
  compare (CTyArr _ _) _ = LT
  compare _ (CTyArr _ _) = GT
  compare (CTyAll _ _ _) _ = LT
  compare _ (CTyAll _ _ _) = GT
  compare (CTyRef _) _ = LT
  compare _ (CTyRef _) = GT
  -- compare (CTyTuple _) _ = LT
  -- compare _ (CTyTuple _) = GT
-}

max_stars_repo_name: minoki/polestar

max_issues_repo_name: minoki/polestar

max_forks_repo_name: minoki/polestar