text:
{-# LANGUAGE LambdaCase, BangPatterns #-}
module Stock where

import Data.Random.Normal
import System.Random
import Data.Random
import Diff
import Numeric.SpecFunctions
import Numeric.LinearAlgebra
import Util

ndigits :: Int
ndigits = 4

roundToN :: Int -> Double -> Double
roundToN n number = (fromInteger (round (number * (10^n)))) / (10.0^n)

-- constants
r :: Double
r = 0.06

variance :: Double
variance = 0.2

days :: Int
days = 360

yperiod :: Double
yperiod = (fromIntegral days) / 360.0

tdiff :: Int
tdiff = 40

tstep :: Double
tstep = (fromIntegral tdiff) / 360.0

s0 :: Double
s0 = 36.0000

k :: Double
k = 40.000

d1 :: Double
d1 = (log (s0 / k) + (r + variance / 2) * yperiod) / (sqrt (variance * yperiod))

d2 :: Double
d2 = d1 - sqrt (variance * yperiod)

-- c = s_0 * N(D_1) - exp (-r * T) * K * N(D_2)
n :: Double -> Double
n = cdf StdNormal

c :: Double
c = s0 * (n d1) - (exp (-r * yperiod)) * k * (n d2)

p :: Double
p = c + k * exp (-r * yperiod) - s0

--- basis functions
l :: Int -> String -> DiffLang
l n x
  | n >= 0    = (Exp ((Const (-0.5)) :*: (Var x))) :*: (Exp (Var x)) :*: (Const (1/(factorial n))) :*: (derivN n (laguerre n x) x)
  | otherwise = error "undefined basis function"

rzs :: Integer -> [Double]
rzs rand = normals' (0.0,1.0) (mkStdGen (fromIntegral rand))

sRand :: Integer -> Int -> Double
sRand j n =
  let zs = rzs j
  in s zs n

s :: [Double] -> Int -> R
s zs 0 = s0
s zs n
  | n > 0 = roundToN ndigits $ (s zs (n - 1)) * (exp $ (r - variance/2) * tstep + ((sqrt (variance * tstep)) * (zs!!n)))

cf :: [R] -> [R]
cf = map (\stdj -> max (k - stdj) 0)

cv :: R -> R
cv = (* (exp (-r * tstep)))

eval :: DiffLang -> R -> R
eval (Var _) x = x
eval (Const n) _ = n
eval (e1 :+: e2) x = eval e1 x + eval e2 x
eval (e1 :*: e2) x = eval e1 x * eval e2 x
eval (e1 :^: e2) x = eval e1 x ** eval e2 x
eval (Exp e) x     = exp (eval e x)

-- compute coefficients {beta_i} using (n-1)'th day's cash as predicator values,
-- n'th days cv as ys
regColN :: ([Vector R], [Vector R]) -> Int -> Matrix R
regColN (paths, cfs) n
  | length paths > n && n > 1 = regCol n laguerreBasis s t
  | otherwise = error "too big a N"
    where st :: ([R], [R], [R])
          st = unzip3 $ filter (\(a,b,c) -> a /= 0) $ zip3 (toList (cfs!!(n-1))) (toList (paths!!(n-1))) (toList (cfs!!n))
          s :: Vector R
          s = fromList $ sndOf3 st
          t :: Vector R
          t = fromList $ thdOf3 st

laguerreJ :: Int -> R -> R
laguerreJ j = eval (laguerre j "X")

laguerreBasis :: [R -> R]
laguerreBasis = map laguerreJ [0..(days `div` tdiff)]

regCol :: Int -> [R -> R] -> Vector R -> Vector R -> Matrix R -- n x 1 result
regCol n basis s t = linearSolveLS a b
  where a :: Matrix R
        -- a = fromColumns $ map (\j -> (fromList . map (basis!!j) . toList) s) [0..n-1]
        a = fromRows $ map (\x -> (fromList . map (\j -> (basis!!j) x)) [0..n-1]) (toList s)
        b :: Matrix R
        b = (asColumn . fromList . map cv) (toList t)

evalPoly :: R -> [R] -> R
evalPoly x betas = foldl (\acc beta -> acc + beta * x) 0 betas

evalPolyL :: [R] -> [R] -> [R]
evalPolyL xs betas = map (\x -> evalPoly x betas) xs

epsilon :: Double
epsilon = 0.00000001

cfP :: ([Vector R], [Vector R]) -> Int -> [Vector R]
cfP (paths, cfs) n
  | length paths > n && n > 1 = [newCFNP, finalCFN]
  | otherwise = error "undefined at the date!"
  where cfn = map (\p -> if fst p == 0 then 0 else snd p) $ zip oldCFP predCF
        finalCFN = fromList $ map (\(s,t) -> if s <= epsilon then 0.0 else t) (zip newCFN oldCFN)
        newCFN = map (\(s,t) -> if t >= s then t else 0.0) (zip oldCFP cfn)
        newCFNP = fromList $ map (\(s,t) -> if t >= s then 0.0 else s) (zip oldCFP cfn)
        -- newCFNP = fromList $ map (\(s,t) -> if t <= epsilon then 0.0 else s) (zip oldCFP cfn)
        oldCFN = toList $ cfs!!n
        oldCFP = toList $ cfs!!(n-1)
        predCF = map (\x -> (foldl (\acc -> \case { (n,beta) -> beta * (laguerreBasis!!n) x + acc }) 0 betas)) (toList (paths!!(n-1)))
        betas :: [(Int,R)]
        betas = zip [0..n-1] ((toList . head . toColumns . regColN (paths, cfs)) n)

cfGenRStep :: ([Vector R], [Vector R]) -> Int -> [Vector R]
cfGenRStep (pathsC, cfsC) n = take (n - 1) cfsC ++ cfP (pathsC, cfsC) n ++ drop (n + 1) cfsC

cfGenD :: (Matrix R, Matrix R) -> [Matrix R]
cfGenD (paths, cfs) = map (fromColumns) $ foldl (\cfsCC n ->
                               let cfs' = head cfsCC
                               in (take (n - 1) cfs' ++ cfP (pathsC, cfs') n ++ drop (n + 1) cfs') : cfsCC)
                      [toColumns cfs] (reverse [2..(days `div` tdiff)])
  where pathsC = toColumns paths

cfGenR :: (Matrix R, Matrix R) -> [Vector R]
cfGenR (paths, cfs) = foldl (\cfs' n -> take (n - 1) cfs' ++ cfP (pathsC, cfs') n ++ drop (n + 1) cfs') cfsC (reverse [2..(days `div` tdiff)])
  where pathsC = toColumns paths
        cfsC = toColumns cfs

tableaux :: Integer -> Integer -> [[R]]
tableaux rand size =
  map (\r -> map (sRand r) [0..(days `div` tdiff)]) ss
  where ss :: [Integer]
        ss = seeds rand size

seeds :: Integer -> Integer -> [Integer]
seeds rand length = take (fromIntegral length) $ map (\x -> (ceiling (abs x * 100000000000))) (rzs rand)

tableauxM :: Integer -> Integer -> (Matrix R, Matrix R)
tableauxM rand size = (matrix (days `div` tdiff + 1) paths, matrix (days `div` tdiff + 1) cfs)
  where paths = concat (tableaux rand size)
        cfs   = cf paths

cfGen :: (Matrix R, Matrix R) -> Matrix R
cfGen = fromColumns . tail . cfGenR

numNonZ :: Matrix R -> Int
numNonZ = sum . map (length . (filter (/= 0)) . toList) . toColumns

aOptionValue :: Matrix R -> [R]
aOptionValue m = map (\lst -> if null lst then 0.0 else (exp (-r * (fromIntegral (fst (last lst))) * tstep)) * (snd (last lst)))
                 $ map (filter (\case {(n,x) -> x /= 0}) . zip [1..(days `div` tdiff)] . toList) (toRows m)

eOptionValue :: Matrix R -> [R]
eOptionValue = map (* (exp (-r * (fromIntegral (days `div` tdiff + 1)) * tstep))) . cf . toList . last . toColumns

cStarMin :: Integer -> (Matrix R, Matrix R) -> R
cStarMin p (cf,m) = (/dist) . sum $ map (\i -> (vs!!(i-1) - vbar) * (qs!!(i-1) - qbar))  [1..(fromIntegral p)]
  where vs = aOptionValue cf
        qs = eOptionValue m
        vbar = sum vs / (fromIntegral (length vs))
        qbar = sum qs / (fromIntegral (length qs))
        dist = sum $ map (\i -> (qs!!(i-1) - qbar)**2) [1..(fromIntegral p)]

rbar :: Integer -> (Matrix R, Matrix R) -> R -> R
rbar pass (cf,m) c = (/(fromIntegral pass)) . sum $ map (\i -> (vs!!(i-1)) - c * ((qs!!(i-1)) - p)) [1..(fromIntegral pass)]
  where vs = aOptionValue cf
        qs = eOptionValue m
        qbar = sum qs / (fromIntegral (length qs))

output :: Integer -> Integer -> R
output rand size = rbar size (cfs, fst m) (cStarMin size (cfs, fst m))
  where cfs = cfGen m
        m = tableauxM rand size

genxs :: Integer -> [Double]
genxs n = map (\i -> let x = aOptionValue (cfGen (tableauxM i 100))
                     in sum x / (fromIntegral (length x)) ) [100..(100+n)]

max_stars_repo_name: TxmszLou/lsmc-cv

max_issues_repo_name: TxmszLou/lsmc-cv

max_forks_repo_name: TxmszLou/lsmc-cv


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE ConstrainedClassMethods #-}
{-# LANGUAGE QuantifiedConstraints #-}

module MLP.Network (
   AllCon,
   ArrList(..),
   arrLast,
   Weights(..),
   Net(..),
   scalarR,
   scalarL,
   MLP(..),
   Learn(..),
   NetInstance
) where

import GHC.TypeLits ( KnownNat, Nat, natVal )
import Numeric.LinearAlgebra.Static
    ( R,
      L,
      (#>),
      (<·>),
      matrix,
      uniformSample,
      unrow,
      vector,
      tr,
      Domain(outer) )
import Fcf ( Eval )
import Fcf.Data.List (Snoc, Last)
import Data.Singletons
    ( Proxy(Proxy), SingI(sing), SingKind(fromSing) )
import Data.Kind (Constraint, Type)
import Control.Monad.Random.Class (MonadRandom(..))
import Numeric.Natural ( Natural )
import Unsafe.Coerce ( unsafeCoerce )
import Data.Proxy ( Proxy(Proxy) )

type family AllCon (f :: k -> Constraint) (ks :: [k]) :: Constraint where
   AllCon _1 '[] = () :: Constraint
   AllCon f (k ': ks) = (f k, AllCon f ks) :: Constraint

data ArrList (e :: Nat -> Type) (a :: [Nat]) where
   Nil :: ArrList e '[]
   Cons :: (KnownNat o, Show (e o)) => e o -> ArrList e os -> ArrList e (o ': os)

arrLast :: forall (o :: Nat) (e :: Nat -> Type) (a :: [Nat]). (AllCon KnownNat a, KnownNat o)
           => ArrList e a -> Maybe (e o)
arrLast Nil = Nothing
arrLast (Cons a as@(Cons _ _)) = arrLast as
arrLast (Cons (a :: e o1) Nil)
   | natVal (Proxy @o1) == natVal (Proxy @o) = Just (unsafeCoerce a)
   | otherwise = Nothing

instance Show (ArrList e a) where
   show Nil = "Nil"
   show (Cons a x) = "Cons " ++ show a ++ show x

data Weights (i :: Nat) (o :: Nat) = Weights {
   inputs :: !(L o i),
   bias :: !(R o)
}

instance (KnownNat i, KnownNat  o) => Show (Weights i o) where
   show (Weights i b) = "\n[ Weights \n" ++ show i ++ "\n, Bias \n" ++ show b ++ "]\n"

rndWeights :: (KnownNat i, KnownNat o, MonadRandom m) => m (Weights i o)
rndWeights = do
   (s1,s2) <- (,) <$> getRandom <*> getRandom

   let w = uniformSample s1 (-1) 1
       b = unrow $ uniformSample s2 (-1) 1

   pure (Weights w b)

data Net (i :: Nat) (hs :: [Nat]) (o :: Nat) where
   SLayer :: !(Weights i o) -> Net i '[] o
   MLayers :: !(Weights i h) -> !(Net h hs o) -> Net i (h ': hs) o

instance AllCon KnownNat (i ': o ': hs) => Show (Net i hs o) where
   show (SLayer l) = "SLayer ->" ++ show l ++ "\n"
   show (MLayers l n) = "MLayers ->" ++ show l ++ "\n" ++ show n

scalarR :: forall n. KnownNat n => Double -> R n
scalarR d = vector (replicate numElems d)
   where
      numElems = fromIntegral . natVal $ Proxy @n

scalarL :: forall n m. (KnownNat n, KnownNat m) => Double -> L n m
scalarL d = matrix (replicate numElems d)
   where
      numElems = product . map fromIntegral . fromSing $ sing @'[n, m]

class MLP a where
   type Topo a :: [Nat]
   type Arr a :: Nat -> Type
   type DimI a :: Nat
   type DimO a :: Nat
   type Layer a :: Nat -> Nat -> Type

   layerOut :: (KnownNat i, KnownNat o)
            => Layer a i o -> Arr a i -> Arr a o

   newLayer :: (KnownNat i, KnownNat o)
               => Double -> Arr a i -> Arr a o -> Layer a i o -> Layer a i o

   netError :: Arr a (DimO a) -> Arr a (DimO a) -> Double

class MLP a => Learn a where
   create :: MonadRandom m => m a

   topo :: a -> [Natural]

   netOut :: Arr a (DimI a) -> a -> ArrList (Arr a) (Topo a)

   mkDeltas :: ArrList (Arr a) (Topo a) -> Arr a (DimO a) -> a -> ArrList (Arr a) (Topo a)

   newNet :: Double -> Arr a (DimI a) -> ArrList (Arr a) (Topo a) -> ArrList (Arr a) (Topo a) -> a -> a

instance AllCon KnownNat (i ': o ': hs)
         => MLP (Net i hs o) where

   type Topo  (Net i hs o) = Eval (Snoc hs o)
   type Arr   (Net i hs o) = R
   type DimI  (Net i hs o) = i
   type DimO  (Net i hs o) = o
   type Layer (Net i hs o) = Weights

   layerOut :: (KnownNat w, KnownNat n) => Weights w n -> R w -> R n
   layerOut l i = sigmoid $ (inputs l #> i) + bias l
      where
         sigmoid x = 1 / (1 + exp (-x))

   newLayer :: (KnownNat w, KnownNat n)
               => Double -> R w -> R n -> Weights w n -> Weights w n
   newLayer lr i delta l = Weights (w - errW) (b - errB)
      where
         w = inputs l
         b = bias l
         errW = scalarL lr * (delta `outer` i)
         errB = scalarR lr * delta

   netError :: KnownNat o => R o -> R o -> Double
   netError expOut actOut = 1 <·> ((expOut - actOut) ** 2)

-- this dummy type can be used to help find the correct MLP instance
-- in functions layerOut, newLayer
type NetInstance = Net 0 '[] 0

instance (KnownNat i, KnownNat o)
         => Learn (Net i '[] o) where

   create :: MonadRandom m => m (Net i '[] o)
   create = SLayer <$> rndWeights

   topo :: Net i '[] o -> [Natural]
   topo _ = [fromSing (sing @o)]

   netOut :: R i -> Net i '[] o -> ArrList R '[o]
   netOut i (SLayer l) = Cons (layerOut @NetInstance l i) Nil

   mkDeltas :: ArrList R '[o] -> R o -> Net i '[] o -> ArrList R '[o]
   mkDeltas (Cons lo Nil) expOut (SLayer l) = Cons (dsn * dso) Nil
      where
         dsn = lo * (1 - lo)
         dso = lo - expOut

   newNet :: Double -> R i -> ArrList R '[o] -> ArrList R '[o] -> Net i '[] o -> Net i '[] o
   newNet lr i (Cons delta Nil) _ (SLayer l) = SLayer (newLayer @NetInstance lr i delta l)

instance (AllCon KnownNat (i ': o ': h ': hs),
         Learn (Net h hs o))
         => Learn (Net i (h ': hs) o) where

   create :: forall i h hs o m. (MonadRandom m, KnownNat i, KnownNat h, Learn (Net h hs o))
             => m (Net i (h ': hs) o)
   create = MLayers <$> (rndWeights @i @h) <*> create @(Net h hs o)

   topo :: Net i (h ': hs) o -> [Natural]
   topo (MLayers _ n) = fromSing (sing @h) : topo n

   netOut :: hhs ~ (h ': hs) => R i -> Net i hhs o -> ArrList R (Topo (Net i hhs o))
   netOut i (MLayers l n) = Cons lo (netOut lo n)
      where
         lo = layerOut @NetInstance l i

   mkDeltas :: hhs ~ (h ': hs) =>
               ArrList R (Topo (Net i hhs o)) -> R o -> Net i hhs o -> ArrList R (Topo (Net i hhs o))
   mkDeltas (Cons lo lor) expOut (MLayers l nxt) =
      case nxt of
         SLayer lnxt -> let dd@(Cons d _) = mkDeltas lor expOut nxt
                            dso = tr (inputs lnxt) #> d
                          in Cons (dsn * dso) dd

         MLayers lnxt _ -> let dd@(Cons d _) = mkDeltas lor expOut nxt
                               dso = tr (inputs lnxt) #> d
                             in Cons (dsn * dso) dd
      where
         dsn = lo * (1 - lo)

   newNet :: Double ->
             R i ->
             ArrList R (Topo (Net i (h ': hs) o)) ->
             ArrList R (Topo (Net i (h ': hs) o)) ->
             Net i (h ': hs) o ->
             Net i (h ': hs) o
   newNet lr i (Cons d ds) (Cons lo los) (MLayers l n) =
      MLayers (newLayer @NetInstance lr i d l) (newNet lr lo ds los n)

layer1 :: Weights 2 3
layer1 = Weights (matrix [1..6]) (vector [1..3])

layer2 :: Weights 3 5
layer2 = Weights (matrix [1..15]) (vector [1..5])

network :: Net 2 '[3] 5
network = MLayers layer1 $ SLayer layer2

ext: hs

max_stars_repo_name: ovanr/Typed-MLP

max_issues_repo_name: ovanr/Typed-MLP

max_forks_repo_name: ovanr/Typed-MLP


text:
{-# LANGUAGE FlexibleContexts #-}
module PlotGreensFunctionR2S1 where

import           Data.Array.Repa         as R
import           Data.Complex
import           Data.List               as L
import           FokkerPlanck.MonteCarlo
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array

{-# INLINE reduceContrast #-}
reduceContrast
  :: (R.Source s Double)
  => Int -> Array s DIM3 Double -> Array D DIM3 Double
reduceContrast idx arr =
  let x = L.head . L.drop idx . L.reverse . L.sort . R.toList $ arr
  in R.map
       (\y ->
          if y >= x
            then x
            else y)
       arr

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:lenStr:initStr:numTrailStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      numTrail = read numTrailStr :: Int
      numThread = read numThreadStr :: Int
  arr <-
    solveMonteCarloR2S1
      numThread
      numTrail
      numPoint
      numPoint
      numOrientation
      sigma
      tao
      numPoint
      ""
  let arr' =
        computeS .
        -- reduceContrast 10 .
        R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . rotate3D $
        arr
      folderPath = "output/app/PlotGreensFunctionR2S1"
  createDirectoryIfMissing True folderPath
  plotImageRepa (folderPath </> printf "GreensR2S1.png") . ImageRepa 8 $
    arr'

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 DeriveFunctor #-}
module Main where

import Lib
import Data.Complex
import Data.Fixed (mod')
import Linear.Vector

data ChebFun a = ChebFun {pow2 :: Int, coefficients :: [a] } deriving (Show, Functor)

-- data SampleFun a =
-- interpolate :: Vector a => SampleFun a -> (Circ -> a) -- need to use linear inteperpolation
{-
instance Functor ChebFun where
	fmap f (ChebFun i cs) = ChebFun i (fmap f cs)
-}

{-
zipWith' f dx dy (x : xs) (y : ys) = (f x y) : (zipWith' f dx dy xs ys)
zipWith' f dx _ [] ys = map (f dx) ys
zipWith' f _ dy xs [] = map (flip f dy) xs
-}

unionWith' f (x : xs) (y : ys) = (f x y) : (unionWith' f xs ys)
unionWith' f [] ys = ys
unionWith' f xs [] = xs

-- Linear package makes scalar multiplication a piece of being a functor.
-- I'm not sure I can make ChebFun a Functor? Well, it can be for the purposes of scalar multiply.
-- But usually we need to FFT it to real space before applying pointwise stuff

instance Additive ChebFun where
	zero = ChebFun 0 []
	liftI2 f (ChebFun n xs) (ChebFun n' ys) = ChebFun (min n n') (zipWith f xs ys)
	liftU2 f (ChebFun n xs) (ChebFun n' ys) = ChebFun (max n n') (unionWith' f xs ys)

-- intance Additive ChebFun where
-- possible conventions: incorrect sized lists are error, everything past nil is 0, everything past last element is a repeat of last constant
-- Last convention allows for efficient representation of constant values. A pure nil is represented as 0. case on (a:[])
-- TrigFun is better name
-- ChebFun b a where ChebFun (b ~ Circ, Num a) =>

-- compose
-- Possible optimization
data Consty f a = NonConsty (f a) | Consty a

{-
we can lift pointwise functions?
cheblift

almost applicative, except that the object inside needs to be foureir transformable.

instance Applicative ChebFun where
	Num a b c => (a -> b -> c) -> ChebFun a -> ChebFun b -> ChebFun c

-}
-- chebcompose :: ChebFun Circ -> ChebFun Circ -> ChebFun Circ

{-
instance Num a => Num (ChubFun a) where
	(ChebFun n xs) + (ChebFun n' ys) = ChebFun (max n n') (zipWithDefault 0 (+) xs ys) -- not right. We want to add with default 0 is [] and then re-shorten.  also the min needs to be a max.
	(ChebFun n xs) * (ChebFun n' ys) = chebify $ (funcify x) * (funcify y) -- all pointwise stuff tends to be in the real domain
	abs x = chebify $ abs . (funicfy x)

instance Vector a => Vector (ChebFun a) where
	smul s = fmap (smul s) coefficients
	vadd x y = x + y

-}

newtype Circ = Circ Double
circ :: Double -> Circ -- smart constructor
circ x = Circ $ mod' (abs x) (2 * pi)

-- https://ro-che.info/articles/2015-12-04-fft
split :: [a] -> ([a], [a])
split = foldr f ([], [])
  where
    f a (r1, r2) = (a : r2, r1)

-- interleave (The opposite of split). Not super sure I got the order right.

-- liquid haskell: interleave split = split interleave
interleave (x : xs, ys) = x : (interleave (ys, xs))

-- only takes 2^n sized lists, or we can assume is zero. Split is incorrect possible then. If not an even number it will swap
fourier :: [Complex Double] -> [Complex Double]
fourier [] = []
fourier [a] = [a]
fourier xs =  (zipWith (+) ffte ffto') <> (zipWith (-) ffte ffto')  where
	        (evens, odds) = split xs
	        ffte = fourier evens
	        n = length ffte
	        ffto = fourier odds
	        ffto' = zipWith (*) ffto (twiddle n)

-- we may want to seperate out the real and imag part
data Fourier = Fourier {k0 :: Double, kcos :: [Double] , ksin :: [Double] }

-- could we get liquid haskell to verify this?

-- I can fourier any sequence Vector a => [a] -> [a]
-- Need scalar multiplication via fourier coefficents.
-- need vector add for + and - of
-- I can make ChebFun a functor. Which can be useful, but perhaps suspicious.

-- instance Num a => Vector a a

fourier' :: [Double] -> [Complex Double]
fourier' [] = []
fourier' [a] = [a :+ 0]
fourier' xs =  (zipWith (+) ffte ffto') <> (zipWith (-) ffte ffto')  where
	        (evens, odds) = split xs
	        ffte = fourier' evens
	        n = length ffte
	        ffto = fourier' odds
	        ffto' = zipWith (*) ffto (twiddle n)

twiddle :: Int -> [Complex Double]
twiddle n = map cis $ iterate (+ (pi / n')) 0 where n' = fromInteger $ toInteger n

-- only takes 2^n sized lists, or we can assume is zero. Split is incorrect possible then. If not an even number it will swap
ifourier :: [Complex Double] -> [Complex Double]
ifourier [] = []
ifourier [a] = [a]
ifourier xs =  (zipWith (+) ffte ffto') <> (zipWith (-) ffte ffto')  where
	        (evens, odds) = split xs
	        ffte = ifourier evens
	        n = length ffte
	        ffto = ifourier odds
	        ffto' = zipWith (*) ffto (map conjugate (twiddle n))
{-
{-# Rewrite (ifourier (fourier x) = x #-}
{-# Rewrite (fourier (ifourier x) = x #-}
	Do I need fourier . ifourier = id ? couldn't hurt I guess.
		These laws aren't strictly true, but pretty close. A nightmare for consistency, but will always make accuracy better.

-}
-- his is some bull.
ifft :: [Complex Double] -> [Complex Double]
ifft (x : xs) = fourier (x : (reverse xs))

fourierGood :: [Complex Double] -> Bool
fourierGood _ = True -- implement the envelope test

envelope :: (Ord a, Num a) => [a] -> [a] -- an envelope that replaces the list with the maximum of the list form then on
envelope = (scanr1 max) . (map abs)  -- (\x acc -> max x acc)

-- plateaujtest j env = env ! j >= env ! j2 * (1 - )   where j2 = 1.25*j + 5
-- hmm a constant 0 may be unhappy with this. probably shouldn't use it.
-- The maximum will be the head of the envelope by the way
-- normalize xs = map (/ m) xs where m = maximum xs

--normalize . envelope

-- The stupidest test. WHat if we are accidentally in a zero of the spectrum
-- maybe lsightly less stupid would be to take the maximum of the last 1/4.
-- These really are garbage. We should actually follow the chebfun prescirption.

simpleTest :: (Fractional a, Ord a) => [a] -> Bool
simpleTest xs = if (last absxs) >= (m * (1e-12)) then False else True where
	    absxs = fmap abs xs
	    m = maximum absxs

-- If we consider iterated ChebFuns, we may want different complexity depending on whether we are at low rank pieces or not.
-- so part of our chebify should recurse inside and chebify for every element.
-- I guess i roughly expect the needed coefficients to form a polytope ?
-- It will be this blob on nenzero coefficients. In higher dimension, surfacc is large, so cutting down on that is helpful.
-- Hmm. This isn't really hierarchical. Maybe some even odd splitting combined with left right splitting in real space?
-- I get a faint sensation I might want to look at a maxi-min? Maybe not.
{-
The new expansion test should decide on both whether to expand internally and at this level.
-- Chebify (Chebify a) => Chebify (Double -> a) where
  chebify f = (fmap chebify) f -

  Chebify Double where
  	chebify = id

-}
-- G(x,x') considered as a 2d function
-- 1/r^2 also should seperate out into peices

-- integrate -- from 0 to x
-- differentiate -- differentiation with resepct to the coordinate
-- functional_differenation -- = differntation with resepct to the value held.
-- sample :: Double -> ChebFun a -> a -- function application
-- fmap diff (For a composed chebfun this idffierentiates the third coordinate)
-- diff2 = diff . diff
-- 3d laplacian.
-- diff2 + fmap diff2 + fmap . fmap diff2

{-
increasingFourier :: [Complex Double] -> (Double -> Double) -> [Complex Double]
increasingFourier evens f | fourierGood evens = increasingFourier evens' f
                          | otherwise = evens where
                                               evens' = (zipWith (+) evens ffto') <> (zipWith (-) evens ffto')
                                               odds = sample f oddpositions
                                               ffto =  fourier odds
                                               twiddle =
                                               ffto' = zipWith      -- if evens acceptable, stop, otherwise

--- increasing fourier is chebify basically
chebify = increasingFourier []
-- sinc interpolation I guess? just take nearest point of reverse fft? Linear interpolation of fft? All of these corresopnd to filters, the choice of
	-- which should not matter much

-- sum_x0  sinc(x-x0) * f(x0)
funcify = \x -> Chebfun
-}

main :: IO ()
main = someFunc


text:
module Test.HUnit.Plus
(
    assertMaybeDouble
  , assertOnFunction
  , assertVector
  , assertMatrix
)
where

import Control.Monad
import Test.HUnit
import Test.HUnit.Approx
import Numeric.LinearAlgebra

assertMaybeDouble :: Maybe Double -> Maybe Double -> Double -> Assertion
assertMaybeDouble Nothing Nothing _ = assertString ""
assertMaybeDouble expected Nothing _ = assertString msg
  where msg = "expected: " ++ show expected ++ "\nbut got: Nothing"
assertMaybeDouble Nothing actual _ = assertString msg
  where msg = "expected: Nothing\nbit got: " ++ show actual
assertMaybeDouble (Just expected) (Just actual) eps = assertApproxEqual "Maybe Double" eps expected actual

assertOnFunction :: (Eq b, Show b) => (a -> b) -> a -> a -> Assertion
assertOnFunction func expected actual = func expected @=? func actual

assertVector :: String -> R -> Vector R -> Vector R -> Assertion
assertVector message eps expected actual =
  let diff = norm_2 (expected-actual)
      msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got" ++ show actual
  in unless (diff < eps) (assertFailure msg)

assertMatrix :: String -> R -> Matrix R -> Matrix R -> Assertion
assertMatrix message eps expected actual =
  let diff = norm_2 (expected-actual)
      msg = message ++ "\nexpected: " ++ show expected ++ "\nbut got" ++ show actual
  in unless (diff < eps) (assertFailure msg)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
-- |
-- Module     : Data.Monoid.Statistics
-- Copyright  : Copyright (c) 2010, Alexey Khudyakov <alexey.skladnoy@gmail.com>
-- License    : BSD3
-- Maintainer : Alexey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability  : experimental
--
module Statistics.Monoid (
    module Statistics.Monoid.Class
  , module Statistics.Monoid.Numeric
  ) where

import Statistics.Monoid.Class
import Statistics.Monoid.Numeric


text:
module LatticeSymmetries.Group
  ( Permutation,
    unPermutation,
    mkPermutation,
    identityPermutation,
    fromGenerators,
    pgLength,
    Symmetry,
    mkSymmetry,
    SymmetriesHeader (..),
    mkSymmetries,
    symmetriesFromHeader,
    Symmetries (..),
    nullSymmetries,
    emptySymmetries,
    symmetriesGetNumberBits,
    symmetriesFromSpec,
  )
where

import Data.Aeson (FromJSON (..), ToJSON (..), object, withObject, (.:), (.=))
import Data.Complex
import qualified Data.List
import Data.Ratio
import qualified Data.Set as Set
import qualified Data.Vector as B
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Storable as S
import Foreign.C.Types (CDouble)
import Foreign.ForeignPtr
import Foreign.Ptr
import Foreign.Storable
import GHC.Exts (IsList (..))
import LatticeSymmetries.Benes
import LatticeSymmetries.Dense
import LatticeSymmetries.FFI
import LatticeSymmetries.IO
import LatticeSymmetries.Utils
import System.IO.Unsafe (unsafePerformIO)
import Prelude hiding (identity, permutations, toList)

newtype PermutationGroup = PermutationGroup (B.Vector Permutation)
  deriving stock (Show, Eq)

getPeriodicity :: Permutation -> Int
getPeriodicity p₀ = go 1 p₀
  where
    identity = identityPermutation (G.length (unPermutation p₀))
    go !n !p
      | p == identity = n
      | otherwise = go (n + 1) (p₀ <> p)

nullPermutationGroup :: PermutationGroup -> Bool
nullPermutationGroup (PermutationGroup gs) = G.null gs

data Symmetry = Symmetry
  { symmetryPermutation :: !Permutation,
    symmetryPhase :: !(Ratio Int)
  }
  deriving stock (Show, Eq, Ord)

symmetrySector :: Symmetry -> Int
symmetrySector s
  | denominator sector == 1 = numerator sector
  | otherwise = error $ "this should not have happened: invalid symmetry: " <> show s
  where
    sector = symmetryPhase s * fromIntegral (getPeriodicity (symmetryPermutation s))

instance FromJSON Symmetry where
  parseJSON = withObject "Symmetry" $ \v -> do
    mkSymmetry <$> v .: "permutation" <*> v .: "sector"

instance ToJSON Symmetry where
  toJSON s = object ["permutation" .= symmetryPermutation s, "sector" .= symmetrySector s]

mkSymmetry :: Permutation -> Int -> Symmetry
mkSymmetry p sector
  | sector < 0 || sector >= periodicity =
    error $
      "invalid sector: " <> show sector <> "; permutation has periodicity " <> show periodicity
  | otherwise = Symmetry p (sector % periodicity)
  where
    periodicity = getPeriodicity p

symmetryNumberSites :: Symmetry -> Int
symmetryNumberSites (Symmetry p _) = G.length (unPermutation p)

modOne :: Integral a => Ratio a -> Ratio a
modOne x
  | x >= 1 = x - fromIntegral (numerator x `div` denominator x)
  | otherwise = x

instance Semigroup Symmetry where
  (<>) (Symmetry pa λa) (Symmetry pb λb) = Symmetry (pa <> pb) (modOne (λa + λb))

data SymmetriesHeader = SymmetriesHeader
  { symmHeaderGroup :: !PermutationGroup,
    symmHeaderNetwork :: !BatchedBenesNetwork,
    symmHeaderCharactersReal :: !(S.Vector Double),
    symmHeaderCharactersImag :: !(S.Vector Double)
  }
  deriving stock (Show, Eq)

getNumberBits :: SymmetriesHeader -> Int
getNumberBits (SymmetriesHeader (PermutationGroup gs) _ _ _)
  | G.null gs = 0
  | otherwise = G.length . unPermutation . G.head $ gs

emptySymmetriesHeader :: SymmetriesHeader
emptySymmetriesHeader = SymmetriesHeader (PermutationGroup G.empty) emptyBatchedBenesNetwork G.empty G.empty

data Symmetries = Symmetries {symmetriesHeader :: !SymmetriesHeader, symmetriesContents :: !(ForeignPtr Cpermutation_group)}
  deriving stock (Show)

instance Eq Symmetries where
  (==) a b = symmetriesHeader a == symmetriesHeader b

instance FromJSON Symmetries where
  parseJSON xs = mkSymmetries <$> parseJSON xs

instance ToJSON Symmetries where
  toJSON symmetries = toJSON (toSymmetryList symmetries)

toSymmetryList :: Symmetries -> [Symmetry]
toSymmetryList (Symmetries (SymmetriesHeader (PermutationGroup gs) _ λsRe λsIm) _) =
  Data.List.zipWith3 _toSymmetry (G.toList gs) (G.toList λsRe) (G.toList λsIm)
  where
    _toSymmetry :: Permutation -> Double -> Double -> Symmetry
    _toSymmetry g λRe λIm
      | r ≈ 1 && s' ≈ fromIntegral (round s' :: Int) = mkSymmetry g (round s')
      | otherwise = error $ "failed to reconstruct Symmetry from " <> show (toList g) <> " and λ = " <> show λRe <> " + " <> show λIm <> "j"
      where
        -- φ = sector / periodicity, but care needs to be taken because φ is approximate
        s' = φ * fromIntegral n
        n = getPeriodicity g
        -- polar returns the phase in (-π, π], we add 2π to get a positive value
        (r, _φ) = polar (λRe :+ λIm)
        φ = _φ / (2 * pi) + (if _φ < 0 then 1 else 0)

pgLength :: PermutationGroup -> Int
pgLength (PermutationGroup gs) = G.length gs

fromGenerators :: (Semigroup a, Ord a) => a -> [a] -> [a]
fromGenerators _ [] = []
fromGenerators identity gs = go Set.empty (Set.singleton identity)
  where
    go !interior !boundary
      | Set.null boundary = Set.toAscList $ interior
      | otherwise = go interior' boundary'
      where
        interior' = interior `Set.union` boundary
        boundary' = Set.fromList [h <> g | h <- Set.toList boundary, g <- gs] Set.\\ interior'

nullSymmetries :: Symmetries -> Bool
nullSymmetries = nullPermutationGroup . symmHeaderGroup . symmetriesHeader

emptySymmetries :: Symmetries
emptySymmetries = symmetriesFromHeader emptySymmetriesHeader

symmetriesGetNumberBits :: Symmetries -> Int
symmetriesGetNumberBits = getNumberBits . symmetriesHeader

mkSymmetries :: [Symmetry] -> Symmetries
mkSymmetries gs = case symmetriesFromHeader <$> mkSymmetriesHeader gs of
  Just s -> s
  Nothing -> error "incompatible symmetries"

mkSymmetriesHeader :: [Symmetry] -> Maybe SymmetriesHeader
mkSymmetriesHeader [] = Just emptySymmetriesHeader
mkSymmetriesHeader gs@(g : _)
  | all ((== n) . symmetryNumberSites) gs = case isConsistent of
    True ->
      let permutations = G.fromList $ symmetryPermutation <$> symmetries
          permGroup = PermutationGroup permutations
          benesNetwork = mkBatchedBenesNetwork $ G.map toBenesNetwork permutations
          charactersReal = G.fromList $ (\φ -> cos (-2 * pi * realToFrac φ)) <$> symmetryPhase <$> symmetries
          charactersImag = G.fromList $ (\φ -> sin (-2 * pi * realToFrac φ)) <$> symmetryPhase <$> symmetries
       in Just $ SymmetriesHeader permGroup benesNetwork charactersReal charactersImag
    False -> Nothing
  | otherwise = error "symmetries have different number of sites"
  where
    n = symmetryNumberSites g
    identity = mkSymmetry (identityPermutation n) 0
    symmetries = fromGenerators identity gs
    isConsistent = all id $ do
      s₁ <- symmetries
      s₂ <- symmetries
      let s₃@(Symmetry p₃ λ₃) = s₁ <> s₂
      pure $
        Set.member s₃ set
          && denominator (λ₃ * fromIntegral (getPeriodicity p₃)) == 1
      where
        set = Set.fromList symmetries

symmetriesFromHeader :: SymmetriesHeader -> Symmetries
symmetriesFromHeader x = unsafePerformIO $ do
  fp <- mallocForeignPtr
  let (DenseMatrix numberShifts numberMasks masks) = bbnMasks $ symmHeaderNetwork x
      shifts = bbnShifts $ symmHeaderNetwork x
      numberBits = getNumberBits x
  withForeignPtr fp $ \ptr ->
    S.unsafeWith masks $ \masksPtr ->
      S.unsafeWith shifts $ \shiftsPtr ->
        S.unsafeWith (symmHeaderCharactersReal x) $ \eigvalsRealPtr ->
          S.unsafeWith (symmHeaderCharactersImag x) $ \eigvalsImagPtr ->
            poke ptr $
              Cpermutation_group
                { cpermutation_group_refcount = 0,
                  cpermutation_group_number_bits = fromIntegral numberBits,
                  cpermutation_group_number_shifts = fromIntegral numberShifts,
                  cpermutation_group_number_masks = fromIntegral numberMasks,
                  cpermutation_group_masks = masksPtr,
                  cpermutation_group_shifts = shiftsPtr,
                  cpermutation_group_eigvals_re = (castPtr :: Ptr Double -> Ptr CDouble) eigvalsRealPtr,
                  cpermutation_group_eigvals_im = (castPtr :: Ptr Double -> Ptr CDouble) eigvalsImagPtr
                }
  pure $ Symmetries x fp
{-# NOINLINE symmetriesFromHeader #-}

symmetryFromSpec :: SymmetrySpec -> Symmetry
symmetryFromSpec (SymmetrySpec p k) =
  mkSymmetry (mkPermutation (G.fromList (toList p))) k

symmetriesFromSpec :: [SymmetrySpec] -> Symmetries
symmetriesFromSpec = mkSymmetries . fmap symmetryFromSpec

max_stars_repo_name: twesterhout/lattice-symmetries-haskell

max_issues_repo_name: twesterhout/lattice-symmetries-haskell

max_forks_repo_name: twesterhout/lattice-symmetries-haskell


text:
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

module ML where

import AI.HNN.FF.Network
import Data.Either (rights)
import Numeric.LinearAlgebra
import System.IO
import System.IO.Unsafe
import Text.ParserCombinators.Parsec

newtype KmH = KmH { kmh :: Double }
   deriving (Num, Floating, Fractional, Real, RealFloat, RealFrac, Eq, Ord)

newtype Rad = Rad { rad :: Double }
   deriving (Num, Floating, Fractional, Real, RealFloat, RealFrac, Eq, Ord)

newtype Force = Force { force :: Double }
   deriving (Num, Floating, Fractional, Real, RealFloat, RealFrac, Eq, Ord)

newtype Sec = Sec { sec :: Double }
   deriving (Num, Floating, Fractional, Real, RealFloat, RealFrac, Eq, Ord)

instance Read KmH where
   readsPrec i = map (\(x,s') -> (KmH x,s')) . readsPrec i

instance Show KmH where
   showsPrec i (KmH d) = showsPrec i d
   show (KmH d) = show d

instance Read Rad where
   readsPrec i = map (\(x,s') -> (Rad x,s')) . readsPrec i

instance Show Rad where
   showsPrec i (Rad d) = showsPrec i d
   show (Rad d) = show d

instance Read Force where
   readsPrec i = map (\(x,s') -> (Force x,s')) . readsPrec i

instance Show Force where
   showsPrec i (Force d) = showsPrec i d
   show (Force d) = show d

instance Read Sec where
   readsPrec i = map (\(x,s') -> (Sec x,s')) . readsPrec i

instance Show Sec where
   showsPrec i (Sec d) = showsPrec i d
   show (Sec d) = show d

data Example = Example { speedX0  :: KmH, speedX1 :: KmH
                       , speedY0  :: KmH, speedY1 :: KmH
                       , angle0   :: Rad, angle1  :: Rad
                       , accel0   :: Force
                       , steer0   :: Force
                       , timeStep :: Sec
                       } deriving Show

deg2rad :: Double -> Rad
deg2rad deg = Rad (deg * 0.017453292519943295)

rad2deg :: Rad -> Double
rad2deg (Rad rad) = rad / 0.017453292519943295

kmh2ms :: KmH -> Double
kmh2ms (KmH kmh) = kmh / 3.6

ms2kmh :: Double -> KmH
ms2kmh ms = KmH (ms * 3.6)

examplesIO :: String -> IO [Example]
examplesIO fn = readFile fn >>= return . map make . rights . map parseLine . lines
   where make :: [Double] -> Example
         make xs = Example { speedX0  = KmH $ xs !! 0, speedX1 = KmH $ xs !! 1
                           , speedY0  = KmH $ xs !! 2, speedY1 = KmH $ xs !! 3
                           , angle0   = Rad $ xs !! 4, angle1  = Rad $ xs !! 5
                           , accel0   = Force $ (xs !! 6) - (xs !! 7)
                           , steer0   = Force $ xs !! 8
                           , timeStep = Sec $ xs !! 9
                           }

         parseLine :: String -> Either ParseError [Double]
         parseLine input = parse lineFmt source input
            where source = "(string: " ++ input ++ ")"

         lineFmt :: GenParser Char st [Double]
         lineFmt = many number

         number :: GenParser Char st Double
         number = do _ <- many space
                     x <- many1 (noneOf "() \n")
                     _ <- many space
                     return (read x)

examples :: [Example]
examples = unsafePerformIO $ examplesIO "bla2"

example2sample :: Example -> Sample Double
example2sample ex = (fromList [{-kmh2ms (speedX0 ex),-} force (accel0 ex)], fromList [kmh2ms (speedX1 ex - speedX0 ex)])

samples :: Samples Double
samples = map example2sample examples

test :: IO ()
test = do   n <- createNetwork 1 [2] 1
            --mapM_ (putStrLn . show . output n tanh . fst) samples
            putStrLn "------------------"
            let n' = trainNTimes 100000 0.8 tanh tanh' n samples
            --mapM_ (putStrLn . show . present n) samples
            putStrLn ("Average error: " ++ show (avgErr n samples))
   where present n (inp, outp) = (toList inp, head $ toList obs, head $ toList outp, head $ diff)
            where obs  = output n tanh inp
                  diff = map (uncurry (-)) (zip (toList obs) (toList outp))
         avgErr n samples = sum (map err samples) / fromIntegral (length samples)
            where err (inp, outp) = abs ((head $ toList $ output n tanh inp) - (head $ toList $ outp)) / abs (head $ toList $ output n tanh inp)

max_stars_repo_name: schwering/golog

max_issues_repo_name: schwering/golog

max_forks_repo_name: schwering/golog


text:
{-# LANGUAGE OverloadedStrings
, DataKinds
, BangPatterns
, TupleSections
, TypeFamilies
, FlexibleContexts #-}
module Pdf.Extract.Spacing
where

-- | This module provides functions for inserting inter-word spaces. A
-- simple rule-based insertion of spaces will work for many texts. But
-- in the case of spaced letters for accentuation, like found in older
-- printings, rule-based insertion might fail. So this module provides
-- an artificial neural network, which can be trained with such cases.

import qualified Data.Text as T
import Data.List
import Data.Maybe
import GHC.Generics hiding (S)
import Control.Lens
import qualified Data.Csv as Csv
import Data.Char
import qualified Data.Vector.Storable as V
import Data.Semigroup
import Grenade
import Grenade.Utils.OneHot
import GHC.TypeNats
import Control.Monad.Random
import Numeric.LinearAlgebra ( maxIndex )
import qualified Numeric.LinearAlgebra.Static as SA
import System.IO
import qualified Data.ByteString.Lazy.Char8 as C8

import Pdf.Extract.Glyph
import Pdf.Extract.Linearize
import Pdf.Extract.Precision

-- * Representation of spaces

-- | Every algorithm for inserting inter-word spaces is considered as
-- a classification algorithm: Each 'Glyph' or 'Char' (or whatever) of
-- a line must be classified either as an item with a space after it,
-- a space before it, spaces around it or as not spaced at all. We use
-- a parametric algebraic data type, `LeftSpacing`, for representing
-- these spacing classes by simply wrapping a glyph or character
-- according to one of these spacing classes.
--
-- We also use 'LeftSpacing' for representing spaces in training and
-- validation data. We can use plain text files with manually verified
-- spaces as training data and transform this plain text into a
-- representation of inter-word spacing using `LeftSpacing`. There are
-- no spaces from the input left in the output characters that are
-- wrapped into 'LeftSpacing'. See 'representSpacesAfter' for such a
-- transformation.
--
-- 'LeftSpacing' is an instance of the functor type class. This
-- enables us to simply use standard 'fmap' to access all the glyph
-- information after space insertion. So in the end, we are free to
-- stack modularized extraction algorithms in whatever order we want.
--
-- Actually only 'SpaceAfter' and 'NoSpace' is used, the others are
-- there for making a conceptionally complete set.
data LeftSpacing a
  = SpaceBefore a               -- ^ Represents a character or 'Glyph'
                                -- with a space before it.
  | SpaceAfter a                -- ^ Represents a character or 'Glyph'
                                -- with a space after it.
  | SpaceAround a               -- ^ Represents a character or 'Glyph'
                                -- with spaces around it.
  | NoSpace a                   -- ^ Represents a character or 'Glyph'
                                -- with neither space before nor after
                                -- it.
  deriving (Show, Eq)

instance Functor LeftSpacing where
  fmap f (SpaceBefore a) = SpaceBefore (f a)
  fmap f (SpaceAfter a) = SpaceAfter (f a)
  fmap f (SpaceAround a) = SpaceAround (f a)
  fmap f (NoSpace a) = NoSpace (f a)

instance Linearizable a => Linearizable (LeftSpacing a) where
  linearize (NoSpace a) = linearizeWithoutState (const ("","")) a
  linearize (SpaceBefore a) = linearizeWithoutState (const (" ","")) a
  linearize (SpaceAfter a) = linearizeWithoutState (const (""," ")) a
  linearize (SpaceAround a) = linearizeWithoutState (const (" "," ")) a

-- | Unwrap a 'LeftSpacing' item.
--
-- Consider using `fmap` instead.
withoutSpace :: LeftSpacing a -> a
withoutSpace (SpaceBefore a) = a
withoutSpace (SpaceAfter a) = a
withoutSpace (NoSpace a) = a

-- | Linearize left spacing text to text.
leftSpacingToText :: [LeftSpacing T.Text] -> T.Text
leftSpacingToText [] = ""
leftSpacingToText ((NoSpace a):xs) = a <> (leftSpacingToText xs)
leftSpacingToText ((SpaceAfter a):xs) = a <> " " <> (leftSpacingToText xs)
leftSpacingToText ((SpaceBefore a):xs) = " " <> a <> (leftSpacingToText xs)

-- | Linearize a line of 'Glyph's to 'T.Text' after space insertion.
leftSpacingGlyphsToText :: Glyph g =>
                           T.Text -- ^ missing text value in glyph data
                        -> [LeftSpacing g]
                        -> T.Text
leftSpacingGlyphsToText missing =
  leftSpacingToText . (map (fmap (fromMaybe missing . text)))

-- * Rule-based insertion of inter-word spaces

-- | Insert spaces on the basis of the size of a glyph. If the
-- distance to the next glyph exceeds the product of the size and a
-- spacing factor, insert a space. In Antiqua scripts (english print)
-- a space was traditionally (at least) a third of the size of a
-- letter. So 1/3 is a good factor to start with.
sizeSpacingFactor :: Glyph g => Double -> [g] -> [LeftSpacing g]
sizeSpacingFactor _ [] = []
sizeSpacingFactor _ (g:[]) = [NoSpace g]
sizeSpacingFactor spacing (g:gn:gs)
  | dist > spacing * (size g)
  = (SpaceAfter g) : sizeSpacingFactor spacing (gn:gs)
  | otherwise
  = (NoSpace g) : sizeSpacingFactor spacing (gn:gs)
  where
    dist = abs((xLeft gn) - (width g) - (xLeft g))

-- | Insert spaces on the basis of the width of a glyph. If the
-- distance to the next glyph exceeds the product of the width and a
-- spacing factor, insert a space.
widthSpacingFactor :: Glyph g => Double -> [g] -> [LeftSpacing g]
widthSpacingFactor _ [] = []
widthSpacingFactor _ (g:[]) = [NoSpace g]
widthSpacingFactor spacing (g:gn:gs)
  | dist > spacing * (width g)
  = (SpaceAfter g) : widthSpacingFactor spacing (gn:gs)
  | otherwise
  = (NoSpace g) : widthSpacingFactor spacing (gn:gs)
  where
    dist = abs((xLeft g) - (xLeft gn))

-- * Artificial neural network for inserting inter-word spaces

-- | Represent training and validation data. This reads a text
-- (`String`) and wraps every character into the `LeftSpacing`
-- functor. Spaces from input string are removed from the output list,
-- but are represented by the functor.
representSpacesAfter :: [Char] -> [LeftSpacing Char]
representSpacesAfter [] = []
representSpacesAfter (x:[])
  | x == ' ' = []
  | otherwise = (NoSpace x):[]
representSpacesAfter (x:x2:xs)
  | x == ' ' = representSpacesAfter (x2:xs) -- drop leading space
  | x2 == ' ' = (SpaceAfter x) : (representSpacesAfter xs)
  | otherwise = (NoSpace x) : (representSpacesAfter (x2:xs))

-- | Helper function for cleaning plain text for training.
cleanForSpaceTraining :: T.Text -> [T.Text]
cleanForSpaceTraining =
  map (T.filter (/=chr 12)) .    -- remove form feeds
  filter (/="") .                -- remove empty lines
  filter (/=(T.pack [chr 12])) . -- remove lines containing form feed
                                 -- only, e.g. last line
  T.lines                        -- split into lines

-- | Validate that the text in a line of input data and in a line of
-- training data is the same.
validateSameText :: Glyph g => [g] -> [LeftSpacing Char] -> Bool
validateSameText glyphs gold =
  (T.concat $ map (fromMaybe "" . text) glyphs) ==
  (T.pack $ map withoutSpace gold)

-- | Generate a data vector from a 'Glyph'.
singleGlyphSpacingVector :: Glyph g => g -> V.Vector Double
singleGlyphSpacingVector g =
  V.fromList $
  [ (width g) / (size g)
  , (size g) / (width g)
  , fromIntegral $ fromEnum $ (size g) < (width g)
  , fromIntegral $ fromEnum $ (size g) == (width g)
  , fromIntegral $ fromEnum $ (size g) > (width g)
  -- , charFeature ord g -- ordinal, causes noise in data and bad results
  , charFeature (fromEnum . isUpper) g
  , charFeature (fromEnum . isLower) g
  , charFeature (fromEnum . (`elem` (",;:.!?"::[Char]))) g
  -- FIXME: make fit for french punctuation
  ]
  where
    charFeature f = fromIntegral . (fromMaybe 0 . (fmap (f . T.head))) . text

-- | Generate data vectors from a line of 'Glyph's by moving a window
-- over the line. The number of preceding and succeeding glyphs,
-- i.e. the size of the window, is set by the first two
-- arguments.
--
-- Missing predecessors (in the beginning of a line) or missing
-- successors (at the end of a line) are represented as missing by
-- one-hot encoded presence/missing vector entries. See
-- 'mkSpacingVector' for details.
mkSpacingVectors
  :: Glyph g =>
     Int                             -- ^ number of preceding glyphs in window
  -> Int                             -- ^ number of succeeding glyphs in window
  -> [g]                             -- ^ the line's glyphs
  -> [V.Vector Double]               -- ^ a vector of float values for each glyph
mkSpacingVectors pre succ gs =
  take (length gs) $ drop pre $
  fst $
  foldr moveWindow ([], (take pre $ repeat Nothing)) $
  ((take pre $ repeat Nothing) ++ (map Just gs) ++ (take succ $ repeat Nothing))
  where
    moveWindow :: Glyph g =>
                  Maybe g
               -> ([V.Vector Double], [Maybe g])
               -> ([V.Vector Double], [Maybe g])
    moveWindow curr (vecs, oldWindow) = ((mkSpacingVector newWindow):vecs, newWindow)
      where
        newWindow = take (pre + succ + 1) (curr:oldWindow)

-- | Make a data vector from a window on some glyphs.
mkSpacingVector :: Glyph g =>
                   [Maybe g]    -- ^ the glyphs in the moving window
                -> V.Vector Double
mkSpacingVector [] = V.empty
mkSpacingVector (g:[]) = V.concat
  [ V.singleton $ fromIntegral $ fromEnum $ isJust g
  , V.singleton $ fromIntegral $ fromEnum $ isNothing g
  , fromMaybe (V.fromList $ take lenSingleVector $ repeat 0.0) $
    fmap singleGlyphSpacingVector g
  ]
  where
    -- SHAPE: adjust to length of singleGlyphSpacingVector
    lenSingleVector = 8
mkSpacingVector (g1:g2:gs) = V.concat
  [ mkSpacingVector (g1:[])
  , (V.singleton $ fromMaybe 0.0 $ (/) <$> dist <*> fmap width g1)
  , mkSpacingVector (g2:gs)
  ]
  where
    dist = spaceBetween <$> g1 <*> g2

-- | Calculate the horizontal space in between two glyphs.
spaceBetween :: (Glyph g) => g -> g -> Double
spaceBetween g1 g2 = abs $ (xLeft g2) - (xRight g1) -- FIXME?

-- | Bring data vectors into shape for training or testing an
-- ANN. This shapes a pair of 'Glyph' vector and 'LeftSpacing'
-- character.
trainingShape :: KnownNat n =>
                 LeftSpacing Char -- ^ label data, i.e. a character wrapped into 'LeftSpacing'
              -> V.Vector Double  -- ^ a data vector
              -> Maybe (S ('D1 n), SpacingOutput)
-- trainingShape (NoSpace _) v = (,) <$> (fromStorable v) <*> oneHot 0
-- trainingShape (SpaceAfter _) v = (,) <$> (fromStorable v) <*> oneHot 1
-- trainingShape (SpaceBefore _) v = (,) <$> (fromStorable v) <*> oneHot 1 -- FIXME: if ever used
-- trainingShape (SpaceAround _) v = (,) <$> (fromStorable v) <*> oneHot 1 -- FIXME: if ever used
trainingShape sp v = (,) <$> (fromStorable v) <*> (oneHot $ leftSpacingToLabel sp)

-- | Make one-hot encoding labels from 'LeftSpacing'.
leftSpacingToLabel :: LeftSpacing a -> Int
leftSpacingToLabel (NoSpace _) = 0
leftSpacingToLabel _ = 1

-- | Return a 'LeftSpacing' constructor for an integer label.
leftSpacingFromLabel :: Int -> (a -> LeftSpacing a)
leftSpacingFromLabel 0 = NoSpace
leftSpacingFromLabel _ = SpaceAfter

-- | Verify training data (or testing data) and generate shaped data
-- for a grenade network using n preceding and m succeeding glyphs,
-- where n and m are given as arguments. This takes a line of text
-- from the set of training data and the corresponding line of glyphs
-- as arguments. The verification tests if both lines contain the same
-- sequence of characters. It is also verified, that a shape was
-- generated for each input glyph.
--
-- If the verification fails, a 'Left' value is returned with an error
-- string.
--
-- Note that this is a generic function. You should use a closure with
-- a type signature like follows to make your shapes of data for a
-- defined number of preceding and succeeding glyphs:
--
-- @myTrainingShapes :: (Glyph g) = T.Text -> [g] -> Either String [(S ('D1 65), SpacingOutput)]@
-- @myTrainingShapes = mkTrainingShapes 0 5@
mkTrainingShapes
  :: (KnownNat n, Glyph g) =>
     Int                    -- ^ number of predecessors in window
  -> Int                    -- ^ number of successors in window
  -> T.Text                 -- ^ a line of text from the training data
  -> [g]                    -- ^ a line of glyphs
  -> Either String [(S ('D1 n), SpacingOutput)]
mkTrainingShapes pre succ txtLine glyphs' =
  join $
  fmap ((ifEitherP
         ((== (length unspacedLine)) . length)
         (const "Error while generating training data: input vector does not fit into shape")
         id) .
        catMaybes .
        (map (uncurry trainingShape)) .
        (zip unspacedLine) .
        (mkSpacingVectors pre succ)) $
  ifEitherP
  ((==(map withoutSpace unspacedLine)) . T.unpack . glyphsTxt)
  (T.unpack .
   T.concat .
   ([ "Error: Line of training data does not match PDF line\n"
       , "Training data:\n"
       , txtLine
       , "\nPDF data:\n"] ++) .
   (:[]) .
   glyphsTxt)
  id
  glyphs
  where
    ifEitherP p l r val
      | p val = Right $ r val
      | otherwise = Left $ l val
    unspacedLine = representSpacesAfter $ T.unpack txtLine
    glyphs = sortOn xLeft glyphs'
    glyphsTxt = T.concat . mapMaybe text -- glyphs

-- | Generate shaped training data with a window showing 2
-- predecessors and 2 successor.
mkTrainingShapes22
  :: (Glyph g) =>
     T.Text                 -- ^ a line of text from the training data
  -> [g]                    -- ^ a line of glyphs
  -> Either String [SpacingRow]
mkTrainingShapes22 = mkTrainingShapes 2 2

-- | Generate shaped data to be run through the network. This is a
-- generic function for generating shaped data with n preceeding and m
-- succeeding glyphs.
mkRunningShapes
  :: (KnownNat n, Glyph g) =>
     Int                    -- ^ number of predecessors in window
  -> Int                    -- ^ number of successors in window
  -> [g]                    -- ^ a line of glyphs
  -> Either String [S ('D1 n)]
mkRunningShapes pre succ glyphs =
  ifEitherP
  ((== (length glyphs)) . length)
  (const "Error while generating data: input vector does not fit into shape")
  id $
  catMaybes $
  map fromStorable $
  mkSpacingVectors pre succ glyphs
  where
    ifEitherP p l r val
      | p val = Right $ r val
      | otherwise = Left $ l val

-- | Generate shaped data for a network with 2 predecessors and 2
-- successors.
mkRunningShapes22 :: Glyph g => [g] -> Either String [SpacingShape]
mkRunningShapes22 = mkRunningShapes 2 2

-- | The shape of data input into the ANN.
--
-- Note that this must be determined from the length of the vector
-- returned by 'singleGlyphSpacingVector', 'mkSpacingVector' and the
-- number of preceding and succeeding glyphs in 'mkTrainingShapes'.
type SpacingShape = S ('D1 54)

-- | The shape of data output from the ANN.
type SpacingOutput = S ('D1 2)

-- | The shape of training data for the ANN.
type SpacingRow = (SpacingShape, SpacingOutput)

-- | The ANN for inserting inter-word spaces.
type SpacingNet
  = Network
    '[ FullyConnected 54 216, Logit,
       FullyConnected 216 216, Logit,
       FullyConnected 216 216, Logit,
       -- FullyConnected 216 216, Logit,
       FullyConnected 216 2, Logit]
    '[ 'D1 54, 'D1 216, 'D1 216,
       'D1 216, 'D1 216,
       'D1 216, 'D1 216,
       -- 'D1 216, 'D1 216,
       'D1 2, 'D1 2]

-- | A net with random weights.
randomSpacingNet :: MonadRandom m => m SpacingNet
randomSpacingNet = randomNetwork

-- | Run a single iteration learning the network.
runSpacingIteration
  :: Handle                     -- ^ a file handle for logging
  -> [SpacingRow]               -- ^ training data
  -> Maybe [SpacingRow]         -- ^ validation data (or Nothing)
  -> LearningParameters         -- ^ learning rate
  -> SpacingNet                 -- ^ the current network
  -> Int                        -- ^ iteration number
  -> IO SpacingNet              -- ^ returns a new network
runSpacingIteration log trainRows validateRows rate net i = do
  let trained' = foldl' (trainEach ( rate { learningRate = learningRate rate * 0.9 ^ i} )) net trainRows
      tRes = fmap (\(rowP,rowL) -> (rowL,) $ runNet trained' rowP) trainRows
      tRes' = fmap (\(S1D label, S1D prediction) -> ( maxIndex (SA.extract label)
                                                    , maxIndex (SA.extract prediction))) tRes
      vRes = fmap (fmap (\(rowP,rowL) -> (rowL,) $ runNet trained' rowP)) validateRows
      vRes' = fmap (fmap (\(S1D label, S1D prediction) ->
                            ( maxIndex (SA.extract label)
                            , maxIndex (SA.extract prediction)))) vRes
  C8.putStr $
    Csv.encode $ (:[]) $
    (precisionLabel 1 tRes'
     & prec_goldStandard .~ (Just "training data")
     & prec_method .~ (Just "2pre2post"))
  case vRes' of
    Nothing -> do
      return ()
    Just res -> do
      C8.putStr $
        Csv.encode $ (:[]) $
        (precisionLabel 1 res
         & prec_goldStandard .~ (Just "validation data")
         & prec_method .~ (Just "2pre2post"))
  return trained'
  where
    trainEach rate' !network (i, o) = train rate' network i o

-- | Run the trained network a line of 'Glyph's. The glyphs should be
-- sorted like in 'linearizeLine'.
runSpacingNetOnLine
  :: (Glyph g) =>
     SpacingNet
  -> [g]
  -> Either String [LeftSpacing g]
runSpacingNetOnLine network glyphs =
  fmap ((map (uncurry (runSpacingNet network))) .
        (zip glyphs)) $
  mkRunningShapes22 glyphs

-- | Run the trained network on input data.
runSpacingNet :: SpacingNet -> a -> SpacingShape -> LeftSpacing a
runSpacingNet net char input =
  spaceFromLabel char $ runNet net input

-- | Evaluate the output from the network.
spaceFromLabel :: a -> SpacingOutput -> LeftSpacing a
spaceFromLabel char (S1D label) = mkSpace char $ maxIndex $ SA.extract label
  where
    mkSpace char 0 = NoSpace char
    mkSpace char 1 = SpaceAfter char

max_stars_repo_name: lueck/scannedb-ok

max_issues_repo_name: lueck/scannedb-ok

max_forks_repo_name: lueck/scannedb-ok


text:
module Corfunc (corfunc,fitdata) where

import Data.Complex
import Numeric.FFT

porod :: Double -> Double -> Double -> Double
porod k sigma q = (k*q**(-4))*exp (-q**2*sigma)

guinier :: Double -> Double -> Double -> Double
guinier a b q = a * exp (b*q**2)

beta :: [Double] -> [Double] -> Double
beta x y = rval x y * std y / std x

alpha :: [Double] -> [Double] -> Double
alpha x y = mean y - beta x y * mean x

mean :: [Double] -> Double
mean xs = sum xs / fromIntegral (length xs)

std :: [Double] -> Double
std xs = sqrt $ mean (map (**2) xs) - mean xs**2

rval :: [Double] -> [Double] -> Double
rval xs ys = sum (zipWith (*) xdiff ydiff) / sqrt (sum (map (**2) xdiff) * sum (map (**2) ydiff))
    where
      diff qs = map (\z -> z - mean qs) qs
      xdiff = diff xs
      ydiff = diff ys

fitguinier :: [Double] -> [Double] -> Double -> Double
fitguinier qs is = guinier a b
    where
      xs = map (**2) qs
      ys = map log is
      a = exp $ alpha xs ys
      b = beta xs ys

fitporod :: [Double] -> [Double] -> Double -> Double
fitporod qs is = porod k sigma
    where
      is' = filter (>0) is
      qs' = filterKey (>0) is qs
      xs = map (**2) qs'
      ys = map log . zipWith (*) is' . map (**4) $ qs'
      k = exp $ alpha xs ys
      sigma = 0 * beta xs ys

filterKey :: (a->Bool) -> [a] -> [b] -> [b]
filterKey f ks = map snd . filter (f.fst) . zip ks

takeKey :: (a->Bool) -> [a] -> [b] -> [b]
takeKey f ks = map snd . takeWhile (f . fst) . zip ks

dropKey :: (a->Bool) -> [a] -> [b] -> [b]
dropKey f ks = map snd . dropWhile (f . fst) . zip ks

fitdata :: [Double] -> [Double] -> Double -> Double
fitdata qs is = let d = interp qs is
                    maxq = 0.04
                    highq = dropWhile (< maxq) qs
                    highi = dropKey (< maxq) qs is
                    fitp = fitporod highq highi
                    minq = 0.0065*3
                    lowq = takeWhile (< minq) qs
                    lowi = takeKey (< minq) qs is
                    fitg = fitguinier lowq lowi
                    s1 = smooth d fitp (maxq,maximum qs)
                    s2 = smooth fitg s1 (minimum qs,minq)
                in
                  s2

-- | Smoothly transition between two functions over a range of functions
smooth :: (Double -> Double) -> (Double -> Double) -> (Double, Double) -> Double -> Double
smooth f g (start, stop) x
    | x <= start = f x
    | x >= stop = g x
    | otherwise = h*g x +(1-h)*f x
    where
      h = 1/(1+(x-stop)**2/(start-x)**2)

interpHead :: ((Double, Double), (Double, Double)) -> Double -> Double
interpHead ((x1,x2),(y1,y2)) q = (x2-q)/(x2-x1)*y1 + (q-x1)/(x2-x1)*y2

interp :: [Double] -> [Double] -> Double -> Double
interp xs ys x =
    let ranges = zip (zip xs (tail xs)) (zip ys (tail ys))
    in
      case dropWhile (\q -> snd (fst q) < x) ranges of
        [] -> last ys
        p:_ -> interpHead p x

dct :: [Double] -> [Double]
dct = map realPart . fft . map (:+ 0)

corfunc :: [Double] -> [Double] -> [(Double,Double)]
corfunc qs is = let fullData = fitdata qs is
                    factor = 16
                    deltaq = (qs!!1-head qs)/factor
                    qgrid = map (*deltaq) [0..factor*2^(floor . log . (/ deltaq) . maximum $ qs)]
                    igrid = zipWith (*) (map (**2) qgrid) $ map fullData qgrid
                    ys = dct igrid
                    xs = map (\q -> 2 * pi * q / deltaq /fromIntegral (length qgrid)) [0..fromIntegral $ length qgrid]
                in
                  zip xs ys

max_stars_repo_name: rprospero/corfunc

max_issues_repo_name: rprospero/corfunc

max_forks_repo_name: rprospero/corfunc


text:
{-# LANGUAGE FlexibleContexts #-}
{- |
Module      : Language.Scheme.Numerical
Copyright   : Justin Ethier
Licence     : MIT (see LICENSE in the distribution)

Maintainer  : github.com/justinethier
Stability   : experimental
Portability : portable

This module implements the numerical tower.
-}

module Language.Scheme.Numerical (
 -- * Generic functions
   numSub
 , numMul
 , numDiv
 , numAdd
 , numMod
 , numRationalize
 , numBoolBinopEq
 , numBoolBinopGt
 , numBoolBinopGte
 , numBoolBinopLt
 , numBoolBinopLte
 , numCast
 , numDenominator
 , numNumerator
 , numInexact2Exact
 , numExact2Inexact
 , num2String
 , unpackNum
 , numericBinop
 -- * Floating point functions
 , numFloor
 , numCeiling
 , numTruncate
 , numRound
 , numExpt
 , numSqrt
 , numExp
 , numLog
 -- * Trigonometric functions
 , numSin
 , numCos
 , numTan
 , numAsin
 , numAcos
 , numAtan
 -- * Complex functions
 , buildComplex
 , numMakePolar
 , numRealPart
 , numImagPart
 , numMagnitude
 , numAngle
 , numMakeRectangular
 -- * Predicates
 , isComplex
 , isReal
 , isRational
 , isInteger
 , isNumber
 , isFloatAnInteger
 , isNumNaN
 , isNumInfinite
 , isNumFinite
 , isNumExact
 , isNumInexact
) where
import Language.Scheme.Types

import Control.Monad.Error
import Data.Char hiding (isNumber)
import Data.Complex
import Data.Fixed
import Data.Ratio
import Numeric
import Text.Printf

-- |A helper function to perform a numeric operation on two values
numericBinop :: (Integer -> Integer -> Integer) -> [LispVal] -> ThrowsError LispVal
numericBinop _ singleVal@[_] = throwError $ NumArgs (Just 2) singleVal
numericBinop op aparams = mapM unpackNum aparams >>= return . Number . foldl1 op

-- - Begin GenUtil - http://repetae.net/computer/haskell/GenUtil.hs
foldlM :: Monad m => (a -> b -> m a) -> a -> [b] -> m a
foldlM f v (x : xs) = (f v x) >>= \ a -> foldlM f a xs
foldlM _ v [] = return v

foldl1M :: Monad m => (a -> a -> m a) -> [a] -> m a
foldl1M f (x : xs) = foldlM f x xs
foldl1M _ _ = error "Unexpected error in foldl1M"
-- end GenUtil

{- FUTURE: as a general comment here, operations need to be more permissive of the
numerical types they accept. Within reason, a user should not have to know
what numerical type they are passing when using these functions -}

-- |Add the given numbers
numAdd :: [LispVal] -> ThrowsError LispVal
numAdd [] = return $ Number 0
numAdd aparams = do
  foldl1M (\ a b -> doAdd =<< (numCast [a, b])) aparams
  where doAdd (List [(Number a), (Number b)]) = return $ Number $ a + b
        doAdd (List [(Float a), (Float b)]) = return $ Float $ a + b
        doAdd (List [(Rational a), (Rational b)]) = return $ Rational $ a + b
        doAdd (List [(Complex a), (Complex b)]) = return $ Complex $ a + b
        doAdd _ = throwError $ Default "Unexpected error in +"

-- |Subtract the given numbers
numSub :: [LispVal] -> ThrowsError LispVal
numSub [] = throwError $ NumArgs (Just 1) []
numSub [Number n] = return $ Number $ -1 * n
numSub [Float n] = return $ Float $ -1 * n
numSub [Rational n] = return $ Rational $ -1 * n
numSub [Complex n] = return $ Complex $ -1 * n
numSub aparams = do
  foldl1M (\ a b -> doSub =<< (numCast [a, b])) aparams
  where doSub (List [(Number a), (Number b)]) = return $ Number $ a - b
        doSub (List [(Float a), (Float b)]) = return $ Float $ a - b
        doSub (List [(Rational a), (Rational b)]) = return $ Rational $ a - b
        doSub (List [(Complex a), (Complex b)]) = return $ Complex $ a - b
        doSub _ = throwError $ Default "Unexpected error in -"

-- |Multiply the given numbers
numMul :: [LispVal] -> ThrowsError LispVal
numMul [] = return $ Number 1
numMul aparams = do
  foldl1M (\ a b -> doMul =<< (numCast [a, b])) aparams
  where doMul (List [(Number a), (Number b)]) = return $ Number $ a * b
        doMul (List [(Float a), (Float b)]) = return $ Float $ a * b
        doMul (List [(Rational a), (Rational b)]) = return $ Rational $ a * b
        doMul (List [(Complex a), (Complex b)]) = return $ Complex $ a * b
        doMul _ = throwError $ Default "Unexpected error in *"

-- |Divide the given numbers
numDiv :: [LispVal] -> ThrowsError LispVal
numDiv [] = throwError $ NumArgs (Just 1) []
numDiv [Number 0] = throwError $ DivideByZero
numDiv [Rational 0] = throwError $ DivideByZero
numDiv [Number n] = return $ Rational $ 1 / (fromInteger n)
numDiv [Float n] = return $ Float $ 1.0 / n
numDiv [Rational n] = return $ Rational $ 1 / n
numDiv [Complex n] = return $ Complex $ 1 / n
numDiv aparams = do
  foldl1M (\ a b -> doDiv =<< (numCast [a, b])) aparams
  where doDiv (List [(Number a), (Number b)])
            | b == 0 = throwError $ DivideByZero
            | (mod a b) == 0 = return $ Number $ div a b
            | otherwise = -- Not an integer
                return $ Rational $ (fromInteger a) / (fromInteger b)
        doDiv (List [(Float a), (Float b)])
            | b == 0.0 = throwError $ DivideByZero
            | otherwise = return $ Float $ a / b
        doDiv (List [(Rational a), (Rational b)])
            | b == 0 = throwError $ DivideByZero
            | otherwise = return $ Rational $ a / b
        doDiv (List [(Complex a), (Complex b)])
            | b == 0 = throwError $ DivideByZero
            | otherwise = return $ Complex $ a / b
        doDiv _ = throwError $ Default "Unexpected error in /"

-- |Take the modulus of the given numbers
numMod :: [LispVal] -> ThrowsError LispVal
numMod [] = return $ Number 1
numMod aparams = do
  foldl1M (\ a b -> doMod =<< (numCast [a, b])) aparams
  where doMod (List [(Number a), (Number b)]) = return $ Number $ mod' a b
        doMod (List [(Float a), (Float b)]) = return $ Float $ mod' a b
        doMod (List [(Rational a), (Rational b)]) = return $ Rational $ mod' a b
        doMod (List [(Complex _), (Complex _)]) = throwError $ Default "modulo not implemented for complex numbers"
        doMod _ = throwError $ Default "Unexpected error in modulo"

-- |Compare a series of numbers using a given numeric comparison
--  function and an array of lisp values
numBoolBinopCompare :: (LispVal
                    -> LispVal -> Either LispError LispVal)
                    -> LispVal -> [LispVal] -> Either LispError LispVal
numBoolBinopCompare cmp n1 (n2 : ns) = do
  (n1', n2') <- numCast' (n1, n2)
  result <- cmp n1' n2'
  case result of
    Bool True -> numBoolBinopCompare cmp n2' ns
    _ -> return $ Bool False
numBoolBinopCompare _ _ _ = return $ Bool True

-- |Numeric equals
numBoolBinopEq :: [LispVal] -> ThrowsError LispVal
numBoolBinopEq [] = throwError $ NumArgs (Just 0) []
numBoolBinopEq (n : ns) = numBoolBinopCompare cmp n ns
  where
    f a b = a == b
    cmp (Number a) (Number b) = return $ Bool $ f a b
    cmp (Float a) (Float b) = return $ Bool $ f a b
    cmp (Rational a) (Rational b) = return $ Bool $ f a b
    cmp (Complex a) (Complex b) = return $ Bool $ f a b
    cmp _ _ = throwError $ Default "Unexpected error in ="

-- |Numeric greater than
numBoolBinopGt :: [LispVal] -> ThrowsError LispVal
numBoolBinopGt [] = throwError $ NumArgs (Just 0) []
numBoolBinopGt (n : ns) = numBoolBinopCompare cmp n ns
  where
    f a b = a > b
    cmp (Number a) (Number b) = return $ Bool $ f a b
    cmp (Float a) (Float b) = return $ Bool $ f a b
    cmp (Rational a) (Rational b) = return $ Bool $ f a b
    cmp _ _ = throwError $ Default "Unexpected error in >"

-- |Numeric greater than equal
numBoolBinopGte :: [LispVal] -> ThrowsError LispVal
numBoolBinopGte [] = throwError $ NumArgs (Just 0) []
numBoolBinopGte (n : ns) = numBoolBinopCompare cmp n ns
  where
    f a b = a >= b
    cmp (Number a) (Number b) = return $ Bool $ f a b
    cmp (Float a) (Float b) = return $ Bool $ f a b
    cmp (Rational a) (Rational b) = return $ Bool $ f a b
    cmp _ _ = throwError $ Default "Unexpected error in >="

-- |Numeric less than
numBoolBinopLt :: [LispVal] -> ThrowsError LispVal
numBoolBinopLt [] = throwError $ NumArgs (Just 0) []
numBoolBinopLt (n : ns) = numBoolBinopCompare cmp n ns
  where
    f a b = a < b
    cmp (Number a) (Number b) = return $ Bool $ f a b
    cmp (Float a) (Float b) = return $ Bool $ f a b
    cmp (Rational a) (Rational b) = return $ Bool $ f a b
    cmp _ _ = throwError $ Default "Unexpected error in <"

-- |Numeric less than equal
numBoolBinopLte :: [LispVal] -> ThrowsError LispVal
numBoolBinopLte [] = throwError $ NumArgs (Just 0) []
numBoolBinopLte (n : ns) = numBoolBinopCompare cmp n ns
  where
    f a b = a <= b
    cmp (Number a) (Number b) = return $ Bool $ f a b
    cmp (Float a) (Float b) = return $ Bool $ f a b
    cmp (Rational a) (Rational b) = return $ Bool $ f a b
    cmp _ _ = throwError $ Default "Unexpected error in <="

-- |Accept two numbers and cast one of them to the appropriate type, if necessary
numCast' :: (LispVal, LispVal) -> ThrowsError (LispVal, LispVal)
numCast' (a@(Number _), b@(Number _)) = return $ (a, b)
numCast' (a@(Float _), b@(Float _)) = return $ (a, b)
numCast' (a@(Rational _), b@(Rational _)) = return $ (a, b)
numCast' (a@(Complex _), b@(Complex _)) = return $ (a, b)
numCast' ((Number a), b@(Float _)) = return $ (Float $ fromInteger a, b)
numCast' ((Number a), b@(Rational _)) = return $ (Rational $ fromInteger a, b)
numCast' ((Number a), b@(Complex _)) = return $ (Complex $ fromInteger a, b)
numCast' (a@(Float _), (Number b)) = return $ (a, Float $ fromInteger b)
numCast' (a@(Float _), (Rational b)) = return $ (a, Float $ fromRational b)
numCast' ((Float a), b@(Complex _)) = return $ (Complex $ a :+ 0, b)
numCast' (a@(Rational _), (Number b)) = return $ (a, Rational $ fromInteger b)
numCast' ((Rational a), b@(Float _)) = return $ (Float $ fromRational a, b)
numCast' ((Rational a), b@(Complex _)) = return $ (Complex $ (fromInteger $ numerator a) / (fromInteger $ denominator a), b)
numCast' (a@(Complex _), (Number b)) = return $ (a, Complex $ fromInteger b)
numCast' (a@(Complex _), (Float b)) = return $ (a, Complex $ b :+ 0)
numCast' (a@(Complex _), (Rational b)) = return $ (a, Complex $ (fromInteger $ numerator b) / (fromInteger $ denominator b))
numCast' (a, b) = case a of
               Number _ -> doThrowError b
               Float _ -> doThrowError b
               Rational _ -> doThrowError b
               Complex _ -> doThrowError b
               _ -> doThrowError a
  where doThrowError num = throwError $ TypeMismatch "number" num

-- |Accept two numbers and cast one of them to the appropriate type, if necessary
numCast :: [LispVal] -> ThrowsError LispVal
numCast [a, b] = do
  (a', b') <- numCast' (a, b)
  pure $ List [a', b']
numCast _ = throwError $ Default "Unexpected error in numCast"

-- |Convert the given number to a rational
numRationalize :: [LispVal] -> ThrowsError LispVal
numRationalize [(Number n)] = return $ Rational $ toRational n
numRationalize [(Float n)] = return $ Rational $ toRational n
numRationalize [n@(Rational _)] = return n
numRationalize [x] = throwError $ TypeMismatch "number" x
numRationalize badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Round the given number
numRound :: [LispVal] -> ThrowsError LispVal
numRound [n@(Number _)] = return n
numRound [(Rational n)] = return $ Number $ round n
numRound [(Float n)] = return $ Float $ fromInteger $ round n
numRound [(Complex n)] = do
  return $ Complex $ (fromInteger $ round $ realPart n) :+ (fromInteger $ round $ imagPart n)
numRound [x] = throwError $ TypeMismatch "number" x
numRound badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Floor the given number
numFloor :: [LispVal] -> ThrowsError LispVal
numFloor [n@(Number _)] = return n
numFloor [(Rational n)] = return $ Number $ floor n
numFloor [(Float n)] = return $ Float $ fromInteger $ floor n
numFloor [(Complex n)] = do
  return $ Complex $ (fromInteger $ floor $ realPart n) :+ (fromInteger $ floor $ imagPart n)
numFloor [x] = throwError $ TypeMismatch "number" x
numFloor badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Take the ceiling of the given number
numCeiling :: [LispVal] -> ThrowsError LispVal
numCeiling [n@(Number _)] = return n
numCeiling [(Rational n)] = return $ Number $ ceiling n
numCeiling [(Float n)] = return $ Float $ fromInteger $ ceiling n
numCeiling [(Complex n)] = do
  return $ Complex $ (fromInteger $ ceiling $ realPart n) :+ (fromInteger $ ceiling $ imagPart n)
numCeiling [x] = throwError $ TypeMismatch "number" x
numCeiling badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Truncate the given number
numTruncate :: [LispVal] -> ThrowsError LispVal
numTruncate [n@(Number _)] = return n
numTruncate [(Rational n)] = return $ Number $ truncate n
numTruncate [(Float n)] = return $ Float $ fromInteger $ truncate n
numTruncate [(Complex n)] = do
  return $ Complex $ (fromInteger $ truncate $ realPart n) :+ (fromInteger $ truncate $ imagPart n)
numTruncate [x] = throwError $ TypeMismatch "number" x
numTruncate badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Sine
numSin :: [LispVal] -> ThrowsError LispVal
numSin [(Number n)] = return $ Float $ sin $ fromInteger n
numSin [(Float n)] = return $ Float $ sin n
numSin [(Rational n)] = return $ Float $ sin $ fromRational n
numSin [(Complex n)] = return $ Complex $ sin n
numSin [x] = throwError $ TypeMismatch "number" x
numSin badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Cosine
numCos :: [LispVal] -> ThrowsError LispVal
numCos [(Number n)] = return $ Float $ cos $ fromInteger n
numCos [(Float n)] = return $ Float $ cos n
numCos [(Rational n)] = return $ Float $ cos $ fromRational n
numCos [(Complex n)] = return $ Complex $ cos n
numCos [x] = throwError $ TypeMismatch "number" x
numCos badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Tangent
numTan :: [LispVal] -> ThrowsError LispVal
numTan [(Number n)] = return $ Float $ tan $ fromInteger n
numTan [(Float n)] = return $ Float $ tan n
numTan [(Rational n)] = return $ Float $ tan $ fromRational n
numTan [(Complex n)] = return $ Complex $ tan n
numTan [x] = throwError $ TypeMismatch "number" x
numTan badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Arcsine
numAsin :: [LispVal] -> ThrowsError LispVal
numAsin [(Number n)] = return $ Float $ asin $ fromInteger n
numAsin [(Float n)] = return $ Float $ asin n
numAsin [(Rational n)] = return $ Float $ asin $ fromRational n
numAsin [(Complex n)] = return $ Complex $ asin n
numAsin [x] = throwError $ TypeMismatch "number" x
numAsin badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Arccosine
numAcos :: [LispVal] -> ThrowsError LispVal
numAcos [(Number n)] = return $ Float $ acos $ fromInteger n
numAcos [(Float n)] = return $ Float $ acos n
numAcos [(Rational n)] = return $ Float $ acos $ fromRational n
numAcos [(Complex n)] = return $ Complex $ acos n
numAcos [x] = throwError $ TypeMismatch "number" x
numAcos badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Arctangent
numAtan :: [LispVal] -> ThrowsError LispVal
numAtan [(Number n)] = return $ Float $ atan $ fromInteger n
numAtan [Number y, Number x] = return $ Float $ phase $ (fromInteger x) :+ (fromInteger y)
numAtan [(Float n)] = return $ Float $ atan n
numAtan [Float y, Float x] = return $ Float $ phase $ x :+ y
numAtan [(Rational n)] = return $ Float $ atan $ fromRational n
numAtan [Rational y, Rational x] = return $ Float $ phase $ (fromRational x) :+ (fromRational y)
numAtan [(Complex n)] = return $ Complex $ atan n
numAtan [x] = throwError $ TypeMismatch "number" x
numAtan badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Take the square root of the given number
numSqrt :: [LispVal] -> ThrowsError LispVal
numSqrt [(Number n)] = if n >= 0 then return $ Float $ sqrt $ fromInteger n
                                 else return $ Complex $ sqrt ((fromInteger n) :+ 0)
numSqrt [(Float n)] = if n >= 0 then return $ Float $ sqrt n
                                else return $ Complex $ sqrt (n :+ 0)
numSqrt [(Rational n)] = numSqrt [Float $ fromRational n]
numSqrt [(Complex n)] = return $ Complex $ sqrt n
numSqrt [x] = throwError $ TypeMismatch "number" x
numSqrt badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Raise the first number to the power of the second
numExpt :: [LispVal] -> ThrowsError LispVal
numExpt [(Number n), (Number p)] = return $ Float $ (fromInteger n) ^ p
numExpt [(Rational n), (Number p)] = return $ Float $ (fromRational n) ^ p
numExpt [(Float n), (Number p)] = return $ Float $ n ^ p
numExpt [(Complex n), (Number p)] = return $ Complex $ n ^ p
numExpt [_, y] = throwError $ TypeMismatch "integer" y
numExpt badArgList = throwError $ NumArgs (Just 2) badArgList

{- numExpt params = do
  foldl1M (\a b -> doExpt =<< (numCast [a, b])) params
  where doExpt (List [(Number a), (Number b)]) = return $ Float $ (fromInteger a) ^ (fromInteger b)
--        doExpt (List [(Rational a), (Rational b)]) = return $ Float $ fromRational $ a ^ b
        doExpt (List [(Float a), (Float b)]) = return $ Float $ a ^ b
--        doExpt (List [(Complex a), (Complex b)]) = return $ Complex $ a ^ b -}

-- |Take the exponent of the given number
numExp :: [LispVal] -> ThrowsError LispVal
numExp [(Number n)] = return $ Float $ exp $ fromInteger n
numExp [(Float n)] = return $ Float $ exp n
numExp [(Rational n)] = return $ Float $ exp $ fromRational n
numExp [(Complex n)] = return $ Complex $ exp n
numExp [x] = throwError $ TypeMismatch "number" x
numExp badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Compute the log of a given number
numLog :: [LispVal] -> ThrowsError LispVal
numLog [(Number n)] = return $ Float $ log $ fromInteger n
numLog [Number n, Number base] = return $ Float $ logBase (fromInteger base) (fromInteger n)
numLog [(Float n)] = return $ Float $ log n
numLog [Float n, Number base] = return $ Float $ logBase (fromInteger base) n
numLog [(Rational n)] = return $ Float $ log $ fromRational n
numLog [Rational n, Number base] = return $ Float $ logBase (fromInteger base) (fromRational n)
numLog [(Complex n)] = return $ Complex $ log n
numLog [Complex n, Number base] = return $ Complex $ logBase (fromInteger base) n
numLog [x] = throwError $ TypeMismatch "number" x
numLog badArgList = throwError $ NumArgs (Just 1) badArgList

-- Complex number functions

-- |Create a complex number
buildComplex :: LispVal
             -- ^ Real part
             -> LispVal
             -- ^ Imaginary part
             -> ThrowsError LispVal
             -- ^ Complex number
buildComplex (Number x) (Number y) = return $ Complex $ (fromInteger x) :+ (fromInteger y)
buildComplex (Number x) (Rational y) = return $ Complex $ (fromInteger x) :+ (fromRational y)
buildComplex (Number x) (Float y) = return $ Complex $ (fromInteger x) :+ y
buildComplex (Rational x) (Number y) = return $ Complex $ (fromRational x) :+ (fromInteger y)
buildComplex (Rational x) (Rational y) = return $ Complex $ (fromRational x) :+ (fromRational y)
buildComplex (Rational x) (Float y) = return $ Complex $ (fromRational x) :+ y
buildComplex (Float x) (Number y) = return $ Complex $ x :+ (fromInteger y)
buildComplex (Float x) (Rational y) = return $ Complex $ x :+ (fromRational y)
buildComplex (Float x) (Float y) = return $ Complex $ x :+ y
buildComplex x y = throwError $ TypeMismatch "number" $ List [x, y]

-- |Create a complex number given its real and imaginary parts
numMakeRectangular :: [LispVal] -> ThrowsError LispVal
numMakeRectangular [x, y] = buildComplex x y
numMakeRectangular badArgList = throwError $ NumArgs (Just 2) badArgList

-- |Create a complex number from its magnitude and phase (angle)
numMakePolar :: [LispVal] -> ThrowsError LispVal
numMakePolar [(Float x), (Float y)] = return $ Complex $ mkPolar x y
numMakePolar [(Float _), y] = throwError $ TypeMismatch "real" y
numMakePolar [x, (Float _)] = throwError $ TypeMismatch "real real" $ x
numMakePolar badArgList = throwError $ NumArgs (Just 2) badArgList

-- |The phase of a complex number
numAngle :: [LispVal] -> ThrowsError LispVal
numAngle [(Complex c)] = return $ Float $ phase c
numAngle [x] = throwError $ TypeMismatch "complex number" x
numAngle badArgList = throwError $ NumArgs (Just 1) badArgList

-- |The nonnegative magnitude of a complex number
numMagnitude :: [LispVal] -> ThrowsError LispVal
numMagnitude [(Complex c)] = return $ Float $ magnitude c
numMagnitude [x] = throwError $ TypeMismatch "complex number" x
numMagnitude badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Retrieve real part of a complex number
numRealPart :: [LispVal] -> ThrowsError LispVal
numRealPart [(Complex c)] = return $ Float $ realPart c
numRealPart [n@(Float _)] = return n
numRealPart [n@(Rational _)] = return n
numRealPart [n@(Number _)] = return n
numRealPart [x] = throwError $ TypeMismatch "complex number" x
numRealPart badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Retrieve imaginary part of a complex number
numImagPart :: [LispVal] -> ThrowsError LispVal
numImagPart [(Complex c)] = do
  let n = imagPart c
      f = Float n
  if isFloatAnInteger f
     then return $ Number $ floor n
     else return f
numImagPart [(Float _)] = return $ Number 0
numImagPart [(Rational _)] = return $ Number 0
numImagPart [(Number _)] = return $ Number 0
numImagPart [x] = throwError $ TypeMismatch "complex number" x
numImagPart badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Take the numerator of the given number
numNumerator :: [LispVal] -> ThrowsError LispVal
numNumerator [n@(Number _)] = return n
numNumerator [(Rational r)] = return $ Number $ numerator r
numNumerator [(Float f)] = return $ Float $ fromInteger . numerator . toRational $ f
numNumerator [x] = throwError $ TypeMismatch "rational number" x
numNumerator badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Take the denominator of the given number
numDenominator :: [LispVal] -> ThrowsError LispVal
numDenominator [Number _] = return $ Number 1
numDenominator [(Rational r)] = return $ Number $ denominator r
numDenominator [(Float f)] = return $ Float $ fromInteger $ denominator $ toRational f
numDenominator [x] = throwError $ TypeMismatch "rational number" x
numDenominator badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Convert an exact number to inexact
numExact2Inexact :: [LispVal] -> ThrowsError LispVal
numExact2Inexact [(Number n)] = return $ Float $ fromInteger n
numExact2Inexact [(Rational n)] = return $ Float $ fromRational n
numExact2Inexact [n@(Float _)] = return n
numExact2Inexact [n@(Complex _)] = return n
numExact2Inexact [badType] = throwError $ TypeMismatch "number" badType
numExact2Inexact badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Convert an inexact number to exact
numInexact2Exact :: [LispVal] -> ThrowsError LispVal
numInexact2Exact [n@(Number _)] = return n
numInexact2Exact [n@(Rational _)] = return n
numInexact2Exact [(Float n)] = return $ Number $ round n
numInexact2Exact [c@(Complex _)] = numRound [c]
numInexact2Exact [badType] = throwError $ TypeMismatch "number" badType
numInexact2Exact badArgList = throwError $ NumArgs (Just 1) badArgList

-- |Convert a number to a string; radix is optional, defaults to base 10
num2String :: [LispVal] -> ThrowsError LispVal
num2String [(Number n)] = return $ String $ show n
num2String [(Number n), (Number radix)] = do
  case radix of
    2 -> do -- Nice tip from StackOverflow question #1959715
             return $ String $ showIntAtBase 2 intToDigit n ""
    8 -> return $ String $ printf "%o" n
    10 -> return $ String $ printf "%d" n
    16 -> return $ String $ printf "%x" n
    _ -> throwError $ BadSpecialForm "Invalid radix value" $ Number radix
num2String [n@(Rational _)] = return $ String $ show n
num2String [(Float n)] = return $ String $ show n
num2String [n@(Complex _)] = return $ String $ show n
num2String [x] = throwError $ TypeMismatch "number" x
num2String badArgList = throwError $ NumArgs (Just 1) badArgList

-- | Determine if the given value is not a number
isNumNaN :: [LispVal] -> ThrowsError LispVal
isNumNaN ([Float n]) = return $ Bool $ isNaN n
isNumNaN _ = return $ Bool False

-- | Determine if number is infinite
isNumInfinite :: [LispVal] -> ThrowsError LispVal
isNumInfinite ([Float n]) = return $ Bool $ isInfinite n
isNumInfinite _ = return $ Bool False

-- | Determine if number is not infinite
isNumFinite :: [LispVal] -> ThrowsError LispVal
isNumFinite ([Number _]) = return $ Bool True
isNumFinite ([Float n]) = return $ Bool $ not $ isInfinite n
isNumFinite ([Complex _]) = return $ Bool True
isNumFinite ([Rational _]) = return $ Bool True
isNumFinite _ = return $ Bool False

-- | Determine if number is exact
isNumExact :: [LispVal] -> ThrowsError LispVal
isNumExact ([Number _]) = return $ Bool True
isNumExact ([Float _]) = return $ Bool False
isNumExact ([Complex _]) = return $ Bool False -- TODO: could be either
isNumExact ([Rational _]) = return $ Bool True
isNumExact _ = return $ Bool False

-- | Determine if number is inexact
isNumInexact :: [LispVal] -> ThrowsError LispVal
isNumInexact ([Number _])   = return $ Bool False
isNumInexact ([Float _])    = return $ Bool True
isNumInexact ([Complex _])  = return $ Bool True
isNumInexact ([Rational _]) = return $ Bool False
isNumInexact _ = return $ Bool False

-- |Predicate to determine if given value is a number
isNumber :: [LispVal] -> ThrowsError LispVal
isNumber ([Number _]) = return $ Bool True
isNumber ([Float _]) = return $ Bool True
isNumber ([Complex _]) = return $ Bool True
isNumber ([Rational _]) = return $ Bool True
isNumber _ = return $ Bool False

-- |Predicate to determine if given number is complex.
--  Keep in mind this does not just look at the types
isComplex :: [LispVal] -> ThrowsError LispVal
isComplex ([Complex _]) = return $ Bool True
isComplex ([Number _]) = return $ Bool True
isComplex ([Rational _]) = return $ Bool True
isComplex ([Float _]) = return $ Bool True
isComplex _ = return $ Bool False

-- |Predicate to determine if given number is a real.
--  Keep in mind this does not just look at the types
isReal :: [LispVal] -> ThrowsError LispVal
isReal ([Number _]) = return $ Bool True
isReal ([Rational _]) = return $ Bool True
isReal ([Float _]) = return $ Bool True
isReal ([Complex c]) = do
  imagPt <- numImagPart [(Complex c)]
  isExact <- isNumExact [imagPt]
  isZero <- numBoolBinopEq [imagPt, (Number 0)]
  case (isExact, isZero) of
    (Bool True, Bool True) -> return $ Bool True
    _ -> return $ Bool False
isReal _ = return $ Bool False

-- |Predicate to determine if given number is a rational.
--  Keep in mind this does not just look at the types
isRational :: [LispVal] -> ThrowsError LispVal
isRational ([Number _]) = return $ Bool True
isRational ([Rational _]) = return $ Bool True
isRational ([Float n]) = return $ Bool $ not $ isInfinite n
isRational _ = return $ Bool False

-- |Predicate to determine if given number is an integer.
--  Keep in mind this does not just look at the types;
--  a floating point input value can return true, for example.
isInteger :: [LispVal] -> ThrowsError LispVal
isInteger ([Number _]) = return $ Bool True
isInteger ([Complex n]) = do
  return $ Bool $ (isFloatAnInteger $ Float $ realPart n) && (isFloatAnInteger $ Float $ imagPart n)
isInteger ([Rational n]) = do
    let numer = abs $ numerator n
    let denom = abs $ denominator n
    return $ Bool $ (numer >= denom) && ((mod numer denom) == 0)
isInteger ([n@(Float _)]) = return $ Bool $ isFloatAnInteger n
isInteger _ = return $ Bool False

-- |A utility function to determine if given value is a floating point
--  number representing an whole number (integer).
isFloatAnInteger :: LispVal -> Bool
isFloatAnInteger (Float n) =
    ((floor n) :: Integer) == ((ceiling n) :: Integer)
isFloatAnInteger _ = False

-- - end Numeric operations section ---

-- |Extract an integer from the given value, throwing a type error if
--  the wrong type is passed.
unpackNum :: LispVal -> ThrowsError Integer
unpackNum (Number n) = return n
unpackNum notNum = throwError $ TypeMismatch "number" notNum

max_stars_repo_name: farrellm/husk-scheme

max_issues_repo_name: farrellm/husk-scheme

max_forks_repo_name: farrellm/husk-scheme