text:
--Utility functions
--Copyright Laurence Emms 2018

module Util (loadIntToRegisters,
             addRegisters,
             ifC,
             whileC,
             doWhileC) where

import Data.Complex
import Data.Bits
import Data.List
import Register
import Instruction
import ClassicalCircuit
import QuantumCircuit

loadIntToRegisters :: (Floating a, Show a, Ord a) => Int -> Int -> [Instruction a]
loadIntToRegisters x r = loadIntBit x r 0 []

loadIntBit :: (Floating a, Show a, Ord a) => Int -> Int -> Int -> [Instruction a] -> [Instruction a]
loadIntBit x r b commands
    | b < 32 = if (testBit x b)
               then (loadIntBit x r (b + 1) (commands ++ [(ITrue (Register (r + b)))]))
               else (loadIntBit x r (b + 1) (commands ++ [(IFalse (Register (r + b)))]))
    | otherwise = commands

addRegisters :: (Floating a, Show a, Ord a) => Int -> Int -> Int -> Int -> Int -> [Instruction a]
addRegisters carry temp a b r
    = (IFalse (Register carry)) :
      (IFalse (Register temp)) :
      (fmap (\(x, y, z) -> (CallCircuit adder
                                       [Right (Register carry),
                                        Right (Register temp),
                                        Right (Register x),
                                        Right (Register y),
                                        Right (Register z)]))
                                       (zip3 [a..(a + 31)] [b..(b + 31)] [r..(r + 31)]))

ifC :: (Floating a, Show a, Ord a) => String -> String -> Classical a -> [Instruction a] -> [Instruction a] -> [Instruction a]
ifC thenLabel endLabel c x y = (JumpWhen thenLabel c) :
                               y ++
                               [(Jump endLabel),
                                (Label thenLabel)] ++
                               x ++
                               [(Label endLabel)]

whileC :: (Floating a, Show a, Ord a) => String -> String -> Classical a -> [Instruction a] -> [Instruction a]
whileC conditionLabel loopLabel c x = [(Jump conditionLabel),
                                       (Label loopLabel)] ++
                                      x ++
                                      [(Label conditionLabel),
                                       (JumpWhen loopLabel c)]

doWhileC :: (Floating a, Show a, Ord a) => String -> Classical a -> [Instruction a] -> [Instruction a]
doWhileC whileLabel c x = (Label whileLabel) :
                          x ++
                          [(JumpWhen whileLabel c)]

max_stars_repo_name: WhatTheFunctional/Hasquil

max_issues_repo_name: WhatTheFunctional/Hasquil

max_forks_repo_name: WhatTheFunctional/Hasquil


text:
{-# LANGUAGE TypeApplications #-}
module ArrayFire.UtilSpec where

import qualified ArrayFire        as A

import           Data.Complex
import           Data.Int
import           Data.Proxy
import           Data.Word
import           Foreign.C.Types
import           System.Directory
import           Test.Hspec

spec :: Spec
spec =
  describe "Util spec" $ do
    it "Should get size of" $ do
      A.getSizeOf (Proxy @ Int) `shouldBe` 8
      A.getSizeOf (Proxy @ Int64) `shouldBe` 8
      A.getSizeOf (Proxy @ Int32) `shouldBe` 4
      A.getSizeOf (Proxy @ Int16) `shouldBe` 2
      A.getSizeOf (Proxy @ Word) `shouldBe` 8
      A.getSizeOf (Proxy @ Word64) `shouldBe` 8
      A.getSizeOf (Proxy @ Word32) `shouldBe` 4
      A.getSizeOf (Proxy @ Word16) `shouldBe` 2
      A.getSizeOf (Proxy @ Word8) `shouldBe` 1
      A.getSizeOf (Proxy @ CBool) `shouldBe` 1
      A.getSizeOf (Proxy @ Double) `shouldBe` 8
      A.getSizeOf (Proxy @ Float) `shouldBe` 4
      A.getSizeOf (Proxy @ (Complex Float)) `shouldBe` 8
      A.getSizeOf (Proxy @ (Complex Double)) `shouldBe` 16
    it "Should get version" $ do
      x <- A.getVersion
      x `shouldBe` (3,6,4)
    it "Should get revision" $ do
      x <- A.getRevision
      x `shouldSatisfy` (not . null)
    it "Should save / read array" $ do
      let arr = A.constant @Int [1,1,1,1] 10
      idx <- A.saveArray "key" arr "file.array" False
      doesFileExist "file.array" `shouldReturn` True
      (`shouldBe` idx) =<< A.readArrayKeyCheck "file.array" "key"
      (`shouldBe` arr) =<< A.readArrayIndex "file.array" idx
      (`shouldBe` arr) =<< A.readArrayKey "file.array" "key"
      removeFile "file.array"

ext: hs

max_stars_repo_name: danielkroeni/arrayfire-haskell

max_issues_repo_name: danielkroeni/arrayfire-haskell

max_forks_repo_name: danielkroeni/arrayfire-haskell


text:
import Data.Map.Strict (Map)
import Numeric.LinearAlgebra
import LeastSquares

type VarName = String
data LinExpr = Var VarName Int
             | ...
infix 4 :==:
data LinContr = LinExpr :==: LinExpr
data QuadExpr = SumSquares LinExpr
              | ProdQuad Double QuadExpr
              | SumQuad QuadExpr QuadExpr

minimize :: QuadExpr -> [LinContr]
                -> Map VarName (Vector R)

max_stars_repo_name: tridao/leastsquares

max_issues_repo_name: tridao/leastsquares

max_forks_repo_name: tridao/leastsquares


text:
module Brain
(
  NetworkWithSamples (..)
, constructNeuralNetwork
, trainN
, activate
)

where

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

import AI.HNN.FF.Network          ( Network
                                  , Samples
                                  , createNetwork
                                  , trainNTimes
                                  , tanh
                                  , tanh'
                                  , output )
import Numeric.LinearAlgebra.Data ( Vector )
import Debug.Trace

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

type Network' = Network Double
type Samples' = Samples Double

data NetworkWithSamples = NetworkWithSamples { network :: Network'
                                             , samples :: Samples' }
                          deriving (Show)

constructNeuralNetwork :: Int -> [Int] -> Int -> Samples' -> IO (NetworkWithSamples)
constructNeuralNetwork i h o s = do
    net <- createNetwork i h o
    return NetworkWithSamples { network = net, samples = s }

trainN :: NetworkWithSamples -> Int -> NetworkWithSamples
trainN nn n = nn { network = nn' }
  where nn' = trainNTimes n 0.05 tanh tanh' (network nn) (samples nn)

activate :: NetworkWithSamples -> Vector Double -> Vector Double
activate n s = output (network n) tanh s


text:
module Test.Statistics where

import Prelude
import qualified Statistics as S
import Test.QuickCheck

runListFold :: S.Statistics t a -> [t] -> a
runListFold = S.runFold foldl

prop_count :: [Bool] -> Property
prop_count l = runListFold (S.count (True ==) 0) l === (length . filter (True ==)) l

main :: IO ()
main = do
    quickCheck prop_count

max_stars_repo_name: maciasoft/wc

max_issues_repo_name: maciasoft/wc

max_forks_repo_name: maciasoft/wc


text:
-- Tasty makes it easy to test your code. It is a test framework that can
-- combine many different types of tests into one suite. See its website for
-- help: <http://documentup.com/feuerbach/tasty>.
import Test.Tasty
-- Hspec is one of the providers for Tasty. It provides a nice syntax for
-- writing tests. Its website has more info: <https://hspec.github.io>.
-- import Test.Tasty.Hspec
import Test.Tasty.HUnit

import ExplicitSimplexStream
import Persistence

import Numeric.LinearAlgebra

-- EXAMPLE STREAMS --

stream1 :: Stream Int
stream1 = addSimplex (addVertex (addVertex (addVertex initializeStream 1) 2) 3) (Simplex [1,2])

stream2 :: Stream Int
stream2 = addVertex (addSimplex (addVertex (addVertex initializeStream 1) 2) (Simplex [2,1])) 3

stream3 :: Stream Int
stream3 = addVertex (addVertex (addVertex (addVertex initializeStream 1) 2) 3) 4

stream4 :: Stream Int
stream4 = initializeStream

stream5 :: Stream Int
stream5 = addSimplex stream4 (Simplex [1,2,3,4])

----

main :: IO ()
main = do
  defaultMain (testGroup "SimplexStream tests" [addSingleVertexTest, initializeStreamTest, streamEqualityTest, streamNumVerticesEmptyTest, streamNumVertices4CellTest, streamGetSizeEmptyTest, streamGetSize4CellTest, streamGetSize4VertexTest, streamToOrderedSimplexListFourVerticesTest, streamToOrderedSimplexListThreeVerticesOneEdgeTest, getBoundaryMapTest, getBoundaryMapTest2, getBoundaryMapTest3, getHomologyDimensionTest, persistenceTest, persistenceTest2, persistenceTest3, persistenceTest4])

addSingleVertexTest :: TestTree
addSingleVertexTest = testCase "Testing addition of single vertex"
  (assertEqual "Should return true for search of vertex 5" (True) (isVertexInStream (addVertex (initializeStream :: Stream Int) 5) 5))

initializeStreamTest :: TestTree
initializeStreamTest = testCase "Testing initialization of stream"
  (assertEqual "Should return Simplices []" (Simplices [Simplex []] :: Stream Int) (initializeStream :: Stream Int))

-- Testing stream equality.
-- Order of Simplex in Stream data type _should not matter_
streamEqualityTest :: TestTree
streamEqualityTest = testCase "Testing quality of streams"
  (assertEqual "Should return Simplices []" (True) (stream1 == stream2))

-- Test numVertices
streamNumVerticesEmptyTest :: TestTree
streamNumVerticesEmptyTest = testCase "Testing number of vertices in empty stream"
  (assertEqual "Should return 0" (0) (numVertices stream4))

streamNumVertices4CellTest :: TestTree
streamNumVertices4CellTest = testCase "Testing number of vertices in a 4-cell"
  (assertEqual "Should return 4" (4) (numVertices stream5))

-- Test getSize
streamGetSizeEmptyTest :: TestTree
streamGetSizeEmptyTest = testCase "Testing get size on empty stream."
  (assertEqual "Should return 1 (the null-cell)" (1) (getSize stream4))

streamGetSize4CellTest :: TestTree
streamGetSize4CellTest = testCase "Testing get size on 4-cell stream."
  (assertEqual "Should return 16" (16) (getSize stream5))

streamGetSize4VertexTest :: TestTree
streamGetSize4VertexTest = testCase "Testing get size on stream with 4 vertices."
  (assertEqual "Should return 5 (four vertices + 1 null cell)." (5) (getSize stream3))

-- Test streamToOrderedSimplexList
streamToOrderedSimplexListFourVerticesTest :: TestTree
streamToOrderedSimplexListFourVerticesTest = testCase "Testing streamToOrderedSimplexList on stream with 4 vertices."
  (assertEqual "Should return object with lengths 0 and 1 with four vertices inside the 1 key." (OrderedSimplexList [SimplexListByDegree 0 [(Simplex [])], SimplexListByDegree 1 [(Simplex [1]), (Simplex [2]), (Simplex [3]), (Simplex [4])]]) (streamToOrderedSimplexList stream3))

streamToOrderedSimplexListThreeVerticesOneEdgeTest :: TestTree
streamToOrderedSimplexListThreeVerticesOneEdgeTest = testCase "Testing streamToOrderedSimplexList on stream with 3 vertices and one edge."
  (assertEqual "Should return object with lengths 0, 1, and 2 with three vertices inside the 1 key and 1 edge inside the 2 key." (OrderedSimplexList [SimplexListByDegree 0 [(Simplex [])], SimplexListByDegree 1 [(Simplex [1]), (Simplex [2]), (Simplex [3])], SimplexListByDegree 2 [(Simplex [1,2])]]) (streamToOrderedSimplexList stream1))

-- Test getBoundaryMap
simplexList1 :: SimplexListByDegree Int
simplexList2 :: SimplexListByDegree Int
simplexList1 = SimplexListByDegree 3 [(Simplex [1,2,3]), (Simplex [2,3,4])]
simplexList2 = SimplexListByDegree 2 [(Simplex [1,2]), (Simplex [2,3]), (Simplex [1,3]), (Simplex [2,4]), (Simplex [3,4])]
getBoundaryMapTest :: TestTree
getBoundaryMapTest = testCase "Testing getBoundaryMap."
  (assertEqual "Should return ..." (fromLists [[1,0],[1,1],[-1,0],[0,-1],[0,1]]) (getBoundaryMap simplexList2 simplexList1))

simplexList3 :: SimplexListByDegree Int
simplexList4 :: SimplexListByDegree Int
simplexList3 = SimplexListByDegree 1 [(Simplex [1]), (Simplex [2]), (Simplex [3]), (Simplex [4])]
simplexList4 = SimplexListByDegree 2 [(Simplex [1,2]), (Simplex [1,3]), (Simplex [1,4]), (Simplex [2,3]), (Simplex [3,4])]
getBoundaryMapTest2 :: TestTree
getBoundaryMapTest2 = testCase "Testing getBoundaryMap."
  (assertEqual "Should return ..." (fromLists [[-1,-1,-1,0,0],[1,0,0,-1,0],[0,1,0,1,-1],[0,0,1,0,1]]) (getBoundaryMap simplexList3 simplexList4))

simplexList5 :: SimplexListByDegree Int
simplexList6 :: SimplexListByDegree Int
simplexList5 = SimplexListByDegree 0 [(Simplex [])]
simplexList6 = SimplexListByDegree 1 [(Simplex [1]), (Simplex [2])]
getBoundaryMapTest3 :: TestTree
getBoundaryMapTest3 = testCase "Testing trivial case for getBoundaryMap."
  (assertEqual "Should return ..." (fromLists [[1, 1]]) (getBoundaryMap simplexList5 simplexList6))

-- Test getHomologyDimension
simplexList7 :: SimplexListByDegree Int
simplexList8 :: SimplexListByDegree Int
simplexList9 :: SimplexListByDegree Int
simplexList7 = SimplexListByDegree 0 [(Simplex [])]
simplexList8 = SimplexListByDegree 1 [(Simplex [1]), (Simplex [2]), (Simplex [3])]
simplexList9 = SimplexListByDegree 2 [(Simplex [1,2]), (Simplex [2,3]), (Simplex [1,3])]
map1 :: Matrix Double
map2 :: Matrix Double
map1 = getBoundaryMap simplexList7 simplexList8
map2 = getBoundaryMap simplexList8 simplexList9
getHomologyDimensionTest :: TestTree
getHomologyDimensionTest = testCase "Testing getHomologyDimension function."
  (assertEqual "Should return ..." (1) (getHomologyDimension map2 map1 0))

-- Test persistence
stream6 :: Stream Int
stream6 = addSimplex (addSimplex (addSimplex initializeStream (Simplex [1,2])) (Simplex [2,3])) (Simplex [1,3])
persistenceTest :: TestTree
persistenceTest = testCase "Testing persistence function."
  (assertEqual "Should return ..." (BettiVector [1,1]) (persistence stream6 1))

persistenceTest2 :: TestTree
persistenceTest2 = testCase "Testing persistence function on filled tetrahedron."
  (assertEqual "Should return ..." (BettiVector [1,0,0,0]) (persistence stream5 1))

persistenceTest3 :: TestTree
persistenceTest3 = testCase "Testing persistence function on four vertices."
  (assertEqual "Should return ..." (BettiVector [4]) (persistence stream3 1))

stream7 :: Stream Int
stream7 = subtractSimplex (addSimplex initializeStream (Simplex [1,2,3,4])) (Simplex [1,2,3,4])
persistenceTest4 :: TestTree
persistenceTest4 = testCase "Testing persistence function."
  (assertEqual "Should return ..." (BettiVector [1,0,1]) (persistence stream7 1))


text:
{-# LANGUAGE FlexibleContexts #-}

module Main (main) where

import Control.Monad.IO.Class (liftIO)
import Data.Complex (Complex)
import Data.Foldable (toList)
import Text.PrettyPrint.Mainland
import Text.PrettyPrint.Mainland.Class

import Spiral
import Spiral.Backend.C
import Spiral.Config
import Spiral.Driver
import Spiral.Exp
import Spiral.OpCount
import Spiral.Program
import Spiral.RootOfUnity
import Spiral.SPL
import Spiral.SPL.Run
import Spiral.Util.Uniq

main :: IO ()
main = defaultMain $ \args -> do
    useComplexType <- asksConfig $ testDynFlag UseComplex
    n <- case args of
           [s] -> return (read s)
           _   -> return 4
    let f = formula n
    pprint f
    if useComplexType
      then toProgram "f" f >>= go
      else toProgram "f" (Re f) >>= go
  where
    go :: (Typed a, Num (Exp a)) => Program a -> Spiral ()
    go prog = do
      pprint prog
      ops <- countProgramOps prog
      resetUnique
      defs <- evalCg $ cgProgram prog
      outp <- asksConfig output
      case outp of
        Nothing -> return ()
        Just{}  -> writeOutput (toList defs)
      liftIO $ putDocLn $
          text "Multiplications:" <+> ppr (mulOps ops) </>
          text "      Additions:" <+> ppr (addOps ops) </>
          text "          Total:" <+> ppr (allOps ops)

-- The SPL formula for which we generate code and count operations.
formula :: Int -> SPL (Exp (Complex Double))
formula _ =
    Pi (L 8 4) ×
    (I 2 ⊗ (Pi (L 4 2) × (I 2 ⊗ F 2 w2) × diag [1.0, 1.0, 1.0, w4] × (F 2 w2 ⊗ I 2))) ×
    diag [1.0, 1.0, 1.0, 1.0, 1.0, w8, w4, w8^(3 :: Int)] ×
    (F 2 w2 ⊗ I 4)
  where
    w2, w4, w8 :: Exp (Complex Double)
    w2 = omega 2
    w4 = omega 4
    w8 = omega 8

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE CPP                        #-}
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE DeriveGeneric              #-}
{-# language FlexibleInstances          #-}
{-# LANGUAGE GADTs                      #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE OverloadedStrings          #-}
{-# LANGUAGE ScopedTypeVariables        #-}
{-# LANGUAGE UndecidableInstances       #-}
{-# LANGUAGE StandaloneDeriving         #-}
{-# language QuasiQuotes                #-}
{-# LANGUAGE TemplateHaskell            #-}
{-# language TypeFamilies               #-}

module Model (
  Val(..),
  ToVal(..),
  FromVal(..),
  Expr(..),
  ModelImage(..),
  Prim1(..),
  Prim2(..),
  Type(..),
  WorkerProfile(..),
  WorkerProfileMap,
  WorkerName(..),
  WorkerProfileId,
  module Type
) where

import Codec.Picture
import Control.Applicative ((<|>))
import Control.DeepSeq
import Control.Monad (mzero)
import Data.ByteString
import Data.Aeson ((.=),(.:))
import qualified Data.ByteString.Char8 as BS
import qualified Data.ByteString.Lazy.Char8 as BL
import qualified Data.ByteString.Base64 as B64
import Data.Complex
import Data.Monoid ((<>))
import Data.Text
import qualified Data.Text as T
import qualified Data.Text.Lazy as TL
import Data.Text.Encoding
import qualified Data.Aeson as A
import qualified Data.Vector as V
#ifndef __GHCJS__
import Database.Groundhog
#endif
import Database.Groundhog.TH hiding (defaultCodegenConfig)
import Generics.SOP
import GHC.Generics
import GHC.TypeLits
import Text.Read
import Text.ParserCombinators.ReadPrec
import qualified Text.PrettyPrint.HughesPJClass as P
import URI.ByteString
import Web.HttpApiData

import EntityID
import Type
import Utils

class ToVal a where
  toVal :: a -> Val

class FromVal a where
  fromVal :: Val -> a

data WorkerName = WorkerName { unWorkerName :: Text }
  deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance NFData WorkerName

data WorkerProfile = WorkerProfile
  { wpName     :: WorkerName
  , wpFunction :: (Text, Type)
  } deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance NFData WorkerProfile

type WorkerProfileMap = EntityMap WorkerProfile
type WorkerProfileId  = EntityID  WorkerProfile

data Expr a = ELit    a Val
            | EVar    a Text
            | ELambda a Text  (Expr a)
            | ERemote a (WorkerProfileId, WorkerProfile, Text)
            | EApp    a (Expr a)  (Expr a)
            | EPrim1  a Prim1 (Expr a)
            | EPrim2  a Prim2 (Expr a) (Expr a)
            deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance Functor Expr where
  fmap f (ELit a v)  = ELit (f a) v
  fmap f (EVar a n)  = EVar (f a) n
  fmap f (ELambda a n b) = ELambda (f a) n (fmap f b)
  fmap f (ERemote a r) = ERemote (f a) r
  fmap f (EApp a eA eB) = EApp (f a) (fmap f eA) (fmap f eB)
  fmap f (EPrim1 a p e) = EPrim1 (f a) p (fmap f e)
  fmap f (EPrim2 a p eA eB) = EPrim2 (f a) p (fmap f eA) (fmap f eB)

instance A.ToJSON a => A.ToJSON (Expr a)
instance A.FromJSON a => A.FromJSON (Expr a)

instance NFData a => NFData (Expr a)

instance NFData Type

data Prim1 = P1Negate | P1Not
  deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance A.ToJSON Prim1
instance A.FromJSON Prim1
instance NFData Prim1

data Prim2 = P2And | P2Or | P2Sum | P2Map | P2Prod
  deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance A.ToJSON Prim2
instance A.FromJSON Prim2
instance NFData Prim2

instance A.ToJSON Type

instance A.FromJSON Type

data Val = VDouble Double
         | VComplex PrimComplex
         | VText Text
         | VImage ModelImage
         | VLabelProbs [(Text, Double)]
         -- TODO: I'm having trouble getting the recursive values into groundhog
         | VList   [Val]
         | VProbabilityDistribution [(Val,Double)]
         | VVec1   (V.Vector Val)
         | VVec2   (V.Vector (Val, Val))
         | VVec3   (V.Vector (Val, Val, Val))
         | VMat2   (V.Vector (V.Vector Val))
         | VMat2C  (V.Vector (V.Vector Val),
                    V.Vector (V.Vector Val),
                    V.Vector (V.Vector Val))
         | VClosure [(Text,Expr Type)] (Expr Type)
         deriving (Eq, Ord, Show, Read, GHC.Generics.Generic)

instance Generics.SOP.Generic Val

newtype PrimComplex = PComplex { getComplex :: Complex Double }
  deriving (Eq, Show, Read, GHC.Generics.Generic, NFData)

instance Ord PrimComplex where
  compare (PComplex p) (PComplex p')
    | realPart p > realPart p' = GT
    | realPart p < realPart p' = LT
    | otherwise = compare (imagPart p) (imagPart p')

instance A.ToJSON PrimComplex where
  toJSON (PComplex c) = A.object ["real" A..= realPart c
                                 ,"imag" A..= imagPart c]

instance A.FromJSON PrimComplex where
  parseJSON (A.Object o) = do
    r <- o A..: "real"
    i <- o A..: "imag"
    return (PComplex (r :+ i))
  parseJSON _ = mzero

-- From basic-sop (not compiling under ghc8)
-- instance NFData Val where
--   rnf = grnf
instance NFData Val

instance A.ToJSON Val
instance A.FromJSON Val

-- | @ModelImage@ is always a base64 encoded tiff with RGBA8 pixels
newtype ModelImage = ModelImage (Image PixelRGBA8)

imageToBytes :: ModelImage -> BS.ByteString
imageToBytes (ModelImage i) = B64.encode . BL.toStrict . encodeTiff $ i

imageFromBytes :: BS.ByteString -> Either String ModelImage
imageFromBytes bs = fmap (ModelImage . convertRGBA8) $ decodeTiff =<< B64.decode bs

instance Show ModelImage where
  show = BS.unpack . imageToBytes

instance Read ModelImage where
  readsPrec _ = readModelImage

readModelImage :: ReadS ModelImage
readModelImage s = case imageFromBytes (BS.pack $ Prelude.dropWhile (== ' ') s) of
  Right m -> [(m,"")]
  Left e -> error e

instance Eq ModelImage where
  a == b =
    A.toJSON a == A.toJSON b

instance Ord ModelImage where
  compare a b = compare (A.encode a) (A.encode b)

instance NFData ModelImage where
  rnf (ModelImage img) = rnf img

instance A.ToJSON ModelImage where
  toJSON (ModelImage img) =
    A.object ["tag"      A..= ("ModelImage" :: String)
             ,"contents" A..= imgString]
    where imgString = A.String . decodeUtf8
                    . B64.encode . BL.toStrict
                    . encodeTiff $ img

instance A.FromJSON ModelImage where
  parseJSON (A.Object o) = do
    t <- o A..: "tag"
    c <- o A..: "contents"
    case (t,c) of
      ("ModelImage" :: String, imgString) ->
        let d = decodeImage =<< B64.decode (encodeUtf8 imgString)
        in case d of
          Right di  -> return $ ModelImage (convertRGBA8 di)
          Left  e    -> error $ "DECODING FAILURE: " ++ e -- mzero
      _ -> mzero

data Tensor = Tensor
  { tDoubleShape :: [Int]
  , tDoubleElems :: V.Vector Double
  }

instance Model.ToVal Int where
  toVal i = Model.VDouble (realToFrac i)

instance Model.ToVal Double where
  toVal = VDouble

instance Model.ToVal Text where
  toVal = VText

instance Model.FromVal Double where
  fromVal (VDouble t) = t
  fromVal e = error $ "Couldn't cast to double: " ++ show e

instance Model.FromVal Text where
  fromVal (VText t) = t
  fromVal e = error $ "Couldn't cast to text: " ++ show e

instance Model.FromVal (Image PixelRGBA8) where
  fromVal (Model.VImage (Model.ModelImage i)) = i
  fromVal x = error $ "Couldn't cast to image: " ++ show x

#ifndef __GHCJS__
mkPersist ghCodeGen [groundhog|
  - primitive: Val
    converter: showReadConverter
|]
#endif

instance A.ToJSON WorkerProfile where
  toJSON (WorkerProfile (WorkerName n) (f,t)) =
    A.object ["name" .= n
             ,"function" .= f
             ,"type" .= t
             ]

instance A.FromJSON WorkerProfile where
  parseJSON (A.Object o) = do
    n <- o .: "name"
    f <- o .: "function"
    t <- o .: "type"
    return $ WorkerProfile n (f,t)

instance FromHttpApiData WorkerName where
  parseUrlPiece = Right . WorkerName

instance A.ToJSON WorkerName where
  toJSON (WorkerName n) = A.String n

instance A.FromJSON WorkerName where
  parseJSON (A.String n) = return $ WorkerName n
  parseJSON _ = mzero

parseWorkerProfile :: Query -> Either String WorkerProfile
parseWorkerProfile q = do
  nm  <- note "No WorkerProfile name"     (lookup "name" ps)
  fn  <- note "No WorkerProfile function" (lookup "function" ps)
  fty <- note "No WorkerProfile type" (lookup "type" ps) >>= (note "No parse" . readMaybe . BS.unpack)
  return $ WorkerProfile (WorkerName $ decodeUtf8 nm) (decodeUtf8 fn, fty)
  where ps = queryPairs q

#ifndef __GHCJS__
-- TODO custom WorkerName instance to avoid showing/reading constructors
mkPersist ghCodeGen [groundhog|
  - primitive: WorkerName
    converter: showReadConverter
  - entity: WorkerProfile
|]
#endif

max_stars_repo_name: CBMM/CBaaS

max_issues_repo_name: CBMM/CBaaS

max_forks_repo_name: CBMM/CBaaS


text:
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}

-----------------------------------------------------------------------------
{-|
Module      : Math.Tensor.LinearAlgebra.Equations
Description : Linear tensor equations.
Copyright   : (c) Nils Alex, 2020
License     : MIT
Maintainer  : nils.alex@fau.de

Linear tensor equations.
-}
-----------------------------------------------------------------------------
module Math.Tensor.LinearAlgebra.Equations
  ( Equation
  , tensorToEquations
  , equationFromRational
  , equationsToSparseMat
  , equationsToMat
  , tensorsToSparseMat
  , tensorsToMat
  , systemRank
  , Solution
  , fromRref
  , fromRrefRev
  , fromRow
  , fromRowRev
  , applySolution
  , solveTensor
  , solveSystem
  , redefineIndets
  ) where

import Math.Tensor
  ( T
  , removeZerosT
  , toListT
  )
import Math.Tensor.LinearAlgebra.Scalar
  ( Poly (Const, Affine, NotSupported)
  , Lin (Lin)
  , polyMap
  , normalize
  )
import Math.Tensor.LinearAlgebra.Matrix
  ( rref
  , isrref
  , verify
  )

import qualified Numeric.LinearAlgebra.Data as HM
  ( Matrix
  , R
  , Z
  , fromLists
  , toLists
  )
import Numeric.LinearAlgebra (rank)

import Data.Maybe (mapMaybe)
import qualified Data.IntMap.Strict as IM
  ( IntMap
  , foldl'
  , map
  , assocs
  , null
  , findWithDefault
  , lookupMax
  , keys
  , fromList
  , (!)
  , difference
  , intersectionWith
  , mapKeys
  , empty
  )
import Data.List (nub, sort)
import Data.Ratio (numerator, denominator)

-- |A linear equation is a mapping from variable
-- indices to coefficients
type Equation a = IM.IntMap a

-- |Extract linear equations from tensor components.
-- The equations are normalized, sorted, and made unique.
tensorToEquations :: Integral a => T (Poly Rational) -> [Equation a]
tensorToEquations = nub . sort . filter (not . IM.null) . fmap (equationFromRational . normalize . snd) . toListT

-- |Extract linear equation with integral coefficients from polynomial
-- tensor component with rational coefficients.
-- Made made integral by multiplying with the @lcm@ of all denominators.
equationFromRational :: forall a.Integral a => Poly Rational -> Equation a
equationFromRational (Affine x (Lin lin))
    | x == 0 = lin'
    | otherwise = error "affine equation not supported for the moment!"
  where
    fac :: a
    fac = IM.foldl' (\acc v -> lcm (fromIntegral (denominator v)) acc) 1 lin
    lin' = IM.map (\v -> fromIntegral (numerator v) * (fac `div` fromIntegral (denominator v))) lin
equationFromRational (Const c)
    | c == 0    = IM.empty
equationFromRational _ = error "equation can only be extracted from linear scalar!"

-- |Convert list of equations to sparse matrix representation of the
-- linear system.
equationsToSparseMat :: [Equation a] -> [((Int,Int), a)]
equationsToSparseMat xs = concat $ zipWith (\i m -> fmap (\(j,v) -> ((i,j),v)) (IM.assocs m)) [1..] xs

-- |Convert list of equations to dense matrix representation of the
-- linear system.
equationsToMat :: Integral a => [Equation a] -> [[a]]
equationsToMat eqns = mapMaybe (\m -> if IM.null m
                                      then Nothing
                                      else Just $ fmap (\j -> IM.findWithDefault 0 j m) [1..maxVar]) eqns
  where
    maxVar = maximum $ mapMaybe (fmap fst . IM.lookupMax) eqns

-- |Extract sparse matrix representation for the linear system given
-- by a list of existentially quantified tensors with polynomial values.
tensorsToSparseMat :: Integral a => [T (Poly Rational)] -> [((Int,Int), a)]
tensorsToSparseMat = equationsToSparseMat . concatMap tensorToEquations

-- |Extract dense matrix representation for the linear system given
-- by a list of existentially quantified tensors with polynomial values.
tensorsToMat :: Integral a => [T (Poly Rational)] -> [[a]]
tensorsToMat = equationsToMat . concatMap tensorToEquations

-- |Rank of a linear system. Uses dense svd provided by hmatrix.
matRank :: forall a.Integral a => [[a]] -> Int
matRank []  = 0
matRank mat = rank (hmat :: HM.Matrix HM.R)
  where
    hmat = HM.fromLists $ fmap (fmap (fromIntegral @a @HM.R)) mat

-- |Rank of the linear system given by a list of existentially
-- quantified tensors with polynomial values.
systemRank :: [T (Poly Rational)] -> Int
systemRank = matRank . tensorsToMat @Int

-- |The solution to a linear system is represented as a list of
-- substitution rules, stored as @'IM.IntMap' ('Poly' 'Rational')@.
type Solution = IM.IntMap (Poly Rational)

-- |Read substitution rules from reduced row echelon form
-- of a linear system.
fromRref :: HM.Matrix HM.Z -> Solution
fromRref ref = IM.fromList assocs
  where
    rows   = HM.toLists ref
    assocs = mapMaybe fromRow rows

fromRrefRev :: HM.Matrix HM.Z -> Solution
fromRrefRev ref = IM.fromList assocs
  where
    rows   = fmap reverse $ HM.toLists ref
    assocs = mapMaybe fromRowRev rows

-- |Read single substitution rule from single
-- row of reduced row echelon form.
fromRow :: forall a.Integral a => [a] -> Maybe (Int, Poly Rational)
fromRow xs = case assocs of
               []             -> Nothing
               [(i,_)]        -> Just (i, Const 0)
               (i, v):assocs' -> let assocs'' = fmap (\(i',v') -> (i', - fromIntegral @a @Rational v' / fromIntegral @a @Rational v)) assocs'
                                 in Just (i, Affine 0 (Lin (IM.fromList assocs'')))
  where
    assocs = filter ((/=0). snd) $ zip [(1::Int)..] xs

fromRowRev :: forall a.Integral a => [a] -> Maybe (Int, Poly Rational)
fromRowRev xs = case assocs of
                  []             -> Nothing
                  [(i,_)]        -> Just (i, Const 0)
                  (i, v):assocs' -> let assocs'' = fmap (\(i',v') -> (i', - fromIntegral @a @Rational v' / fromIntegral @a @Rational v)) assocs'
                                    in Just (i, Affine 0 (Lin (IM.fromList assocs'')))
  where
    assocs = reverse $ filter ((/=0). snd) $ zip [(1::Int)..] xs

-- |Apply substitution rules to tensor component.
applySolution :: Solution -> Poly Rational -> Poly Rational
applySolution s (Affine x (Lin lin))
    | x == 0 = case p of
                 Affine xFin (Lin linFin) -> if IM.null linFin
                                             then Const xFin
                                             else p
                 _ -> p
    | otherwise = error "affine equations not yet supported"
  where
    s' = IM.intersectionWith (\row v -> if v == 0
                                        then error "value 0 encountered in linear scalar"
                                        else polyMap (v*) row) s lin
    lin' = IM.difference lin s
    p0 = if IM.null lin'
         then Const 0
         else Affine 0 (Lin lin')

    p = IM.foldl' (+) p0 s'

    {-
    lin' = IM.foldlWithKey'
             (\lin' i sub ->
                 case Affine 0 (Lin lin') + sub of
                   Affine 0 (Lin lin'') -> IM.delete i lin''
                   Const 0              -> IM.empty
                   _                    -> error "affine equations not yet supported")
             lin s'
    -}
applySolution _ _ = error "only linear equations supported"

-- |Apply substitution rules to all components of a tensor.
solveTensor :: Solution -> T (Poly Rational) -> T (Poly Rational)
solveTensor sol = removeZerosT . fmap (applySolution sol)

-- |Solve a linear system and apply solution to the tensorial
-- indeterminants.
solveSystem ::
     [T (Poly Rational)] -- ^ Tensorial linear system
  -> [T (Poly Rational)] -- ^ List of indeterminant tensors
  -> [T (Poly Rational)] -- ^ Solved indeterminant tensors
solveSystem system indets
    | wrongSolution = error "Wrong solution found. May be an Int64 overflow."
    | otherwise     = indets'
  where
    mat = HM.fromLists $ tensorsToMat @HM.Z system
    ref = rref mat
    wrongSolution = not (isrref ref && verify mat ref)
    sol = fromRref ref
    indets' = fmap (solveTensor sol) indets

-- |Relabelling of the indeterminants present in a list of tensors.
-- Redefines the labels of @n@ indeterminants as @[1..n]@, preserving
-- the previous order.
redefineIndets :: [T (Poly v)] -> [T (Poly v)]
redefineIndets indets = fmap (fmap (\case
                                       Const c -> Const c
                                       NotSupported -> NotSupported
                                       Affine a (Lin lin) ->
                                         Affine a (Lin (IM.mapKeys (varMap IM.!) lin)))) indets
  where
    comps = snd <$> concatMap toListT indets
    vars = nub $ concat $ mapMaybe (\case
                                       Affine _ (Lin lin) -> Just $ IM.keys lin
                                       _                  -> Nothing) comps
    varMap = IM.fromList $ zip vars [(1::Int)..]

max_stars_repo_name: nilsalex/type-tensor

max_issues_repo_name: nilsalex/type-tensor

max_forks_repo_name: nilsalex/type-tensor


text:
module Y2017.M02.D24.Solution where

-- below imports from cabal package 'statistics'

import Statistics.Distribution
import Statistics.Distribution.Binomial

{--
AHA!

Whenever you think you have to invent something, stop, pause, and think, ...

SOMEBODY ELSE HAS ALREADY DONE IT! USE THE DARK SIDE OF THE FORCE AND TOGETHER
WE CAN RULE THE GALAXY AS FATHER AND SON!

See? All you need is love, a psychopathic father that murders young children
at the Jedi Academy, and 'midichlorians.' Lots of 'midichlorians.'

Where was I?

Ah, yes: 'already invented.'

So, we know from Rosalind.Traits that any descendant from an Aa Bb/Aa Bb
pairing will have a 25% probability of itself being Aa Bb (you can verify
this if you'd like), so the question arises, what is the probability of at
least n offspring at generation x (or Generation X, if you prefer) also being
Aa Bb.

This is the cumulative (or actually co-cumulative) binomial probability
distribution.

So.

Today's Haskell problem: with the help of the above imports, create a binomial
probability distribution at 25%, run a scenario where you are in generation
x, where the xth generation has x^2 offspring, and you are looking for at least
n of the offspring to have the desired characteristic (in the abstract case,
the what the characteristic is doesn't matter, we need what the probably outcome
of having it.
--}

type Trials = Int

aAbB :: Trials -> BinomialDistribution
aAbB = (`binomial` 0.25)

genx :: Int -> BinomialDistribution -> Double
genx n bin = sum (map (probability bin) [n .. bdTrials bin])

{--
What is the result for gen = 2, n = at least 1?

>>> genx 1 (aAbB 4)
0.68359375

Now solve for gen = 7, n = at least 35:

>>> genx 35 (aAbB (2 ^ 7))
0.3006022057910504

*WHEW* That was a lot easier than rewriting probability theory by hand!
--}


text:
module Lib.Practice.P0
    (
        f0
    ) where

import Numeric.LinearAlgebra
import Lib.Utility

f0 = do
    -- 配列の生成
    let x = vector [1..3]
    print x

    -- 加減乗除
    let y = vector [2, 4..6]
    print $ x + y
    print $ x - y
    print $ x * y
    print $ x / y

    -- ブロードキャスト
    print $ x / 2

    -- 行列
    let mA = (2><2) [1..4] :: Matrix Double
    print mA
    print $ size mA
    let mB = (2><2) [3, 0, 0, 6]
    print $ mA + mB
    print $ mA * mB

    -- 行列のブロードキャスト
    print $ mA * 10
    print $ mA * row [10, 20]
    -- ↑MatrixとVectorのelement-wiseな計算用の演算子がない？

    -- 要素へのアクセス
    let mX = (3><2) [51, 55, 14, 19, 0, 4 :: Double]
    print $ mX ! 0
    print $ mX ! 0 ! 1
    mapM_ print $ toRows mX
    let flattened = flatten mX
    print $ map (flattened !) $ [0, 2, 4]
    print $ map (==1) $ toList $ flattened .>. 15
    print $ filter (>15) $ toList flattened


text:
{-|
Module: MachineLearning.NeuralNetwork.WeightInitialization
Description: Weight Initialization
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Various Weight Initialization algorithms.
-}

module MachineLearning.NeuralNetwork.WeightInitialization
(
  nguyen
  , he
)

where

import qualified Numeric.LinearAlgebra as LA
import qualified Control.Monad.Random as RndM
import MachineLearning.Types (Matrix)
import MachineLearning.Random (getRandomRMatrixM)

-- | Weight Initialization Algorithm discussed in Nguyen at al.: https://web.stanford.edu/class/ee373b/nninitialization.pdf
-- Nguyen, Derrick, Widrow, B. Improving the learning speed of 2-layer neural networks by choosing initial values of adaptive weights.
-- In Proc. IJCNN, 1990; 3: 21-26.
nguyen :: RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
nguyen (r, c) = do
  let b = LA.konst 0 (1, r)
      eps = (sqrt 6) / (sqrt . fromIntegral $ r + c)
  w <- getRandomRMatrixM r c (-eps, eps)
  return (b, w)

-- | Weight Initialization Algorithm discussed in He at al.: https://arxiv.org/abs/1502.01852
-- Kaiming He, Xiangyu Zhang, Shaoqing Ren, Jian Sun.
-- Delving Deep into Rectifiers: Surpassing Human-Level Performance on ImageNet Classification.
he :: RndM.RandomGen g => (Int, Int) -> RndM.Rand g (Matrix, Matrix)
he (r, c) = do
  let b = LA.konst 0 (1, r)
      eps = sqrt (2/(fromIntegral $ r + c))
  w <- getRandomRMatrixM r c (-eps, eps)
  return (b, w)

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module Main where

import           Control.Monad.ST               ( runST )
import qualified Data.Vector                   as V
import           Statistics.BBVI.Examples

main :: IO ()
main =
  let xs = runST $ genMixedMem
  in  (putStrLn . unlines . fmap show . take 100 . V.toList) xs
        >> putStrLn (show . snd $ mixedMemFit xs)

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
{-# LANGUAGE CPP #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Properties (tests) where

import Control.Applicative as A ((<$>))
import Criterion.Analysis
import Prelude ()
import Prelude.Compat
import Statistics.Types (Sample)
import Test.Tasty (TestTree, testGroup)
import Test.Tasty.QuickCheck (testProperty)
import Test.QuickCheck
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U

instance (Arbitrary a, U.Unbox a) => Arbitrary (U.Vector a) where
  arbitrary = U.fromList A.<$> arbitrary
  shrink    = map U.fromList . shrink . U.toList

outlier_bucketing :: Double -> Sample -> Bool
outlier_bucketing y ys =
  countOutliers (classifyOutliers xs) <= fromIntegral (G.length xs)
  where xs = U.cons y ys

outlier_bucketing_weighted :: Double -> Sample -> Bool
outlier_bucketing_weighted x xs =
  outlier_bucketing x (xs <> G.replicate (G.length xs * 10) 0)

tests :: TestTree
tests = testGroup "Properties" [
    testProperty "outlier_bucketing" outlier_bucketing
  , testProperty "outlier_bucketing_weighted" outlier_bucketing_weighted
  ]

max_stars_repo_name: rubenpieters/criterion

max_issues_repo_name: rubenpieters/criterion

max_forks_repo_name: rubenpieters/criterion


text:
module FET where

import qualified Data.ByteString.Lazy.Char8             as BS
import           Data.Eq
import           Data.Function
import           Data.Functor
import           Data.Int
import           Data.List
import           Data.Ord
import           Prelude                                (Double, fromIntegral,
                                                         (*), (+), (-), (/),
                                                         (^), otherwise, (&&))
import           Statistics.Distribution                (complCumulative,
                                                         cumulative,
                                                         probability)
import           Statistics.Distribution.ChiSquared
import           Statistics.Distribution.Hypergeometric
import           Text.Show

--the Fisher's exact test can be computed from the hypergeometric distribution
--eg. test whether the presence / absence of a gene is statistically significant
--between two bacterial groups
--           Group 1 | Group 2
--          -------------------
-- has gene |   2    |    7
-- no gene  |   5    |    3
--          -------------------
-- the hypergeometric distribution is calculated in Statistics.Distribution.Hypergeometric
-- hypergeometric = m l k
-- From the package description: The parameters of the distribution describe
-- k elements chosen from a population of l, with m elements of one type,
-- and l-m of the other (all are positive integers).
-- In our case:
-- m is the sum of "has factor" instances, so 9
-- l is the total population size, so 17
-- k is the sum of a particular group, so 7 for Group 1
-- let d = distribution 9 17 7
-- This gives us the distribution, but we need to test the probability
-- that the presence of this gene is statistically different, using the two-tail
-- approach.
-- The one-tail probability in our example of Group 1 having 2 members with the
-- gene can be tested using the cumulative distribution function from
-- Statistics.Distribution, our previously calculated distribution "d", and
-- the number of "successes" (x), which is 2.
-- let oneTail = cumulative d 2 = 0.11703002879473468
-- However, we have no previous expectation for this group to have or not have
-- this gene, so what we really want is the two-tailed test.
-- Please see http://www.pmean.com/11/Fishers.html for an excellent,
-- detailed description.
-- In short, we need to sum all of the probabilities that are smaller than or
-- equal to the probability of the observed table (eg. x = 2 in our example).
-- The twoTailValue is calculated by generating all of the probabilities for
-- our sample, filtering the ones that have a probability too large, and
-- summing the remaining values.
--
-- For a total number of observations greater than 1000, we use the Chi Square
-- test, as the FET returns NaN.
-- ChiSqr statistic  for
--           Group 1 | Group 2
--          -------------------
-- has gene |   2    |    7
-- no gene  |   5    |    3
--          -------------------
-- is [(2*3)-(5*7)]^2 * (2+7+5+3) / (2+7)(5+3)(7+3)(2+5)
-- With 1 degree of freedom, we can calculate the Chi Square distribution
-- with chiSquared 1 from Statistics.Distribution
-- To determine the P value of a Chi Square statistic, we use the
-- complCumulative function and the
-- Chi Square distribution with 1 degree of freedom.
-- For Chi Square, the one-tail is simply half the two-tail
-- Statistics.Test.ChiSquared only return "Significant" or "NotSignificant"
-- and we wanted the actual P-value, hence the function.
-- Likewise there is no Statistics.Test.FET

newtype GroupOneA = GroupOneA { unGroupOneA :: Int } deriving (Eq, Show)
newtype GroupOneB = GroupOneB { unGroupOneB :: Int } deriving (Eq, Show)
newtype GroupTwoA = GroupTwoA { unGroupTwoA :: Int } deriving (Eq, Show)
newtype GroupTwoB = GroupTwoB { unGroupTwoB :: Int } deriving (Eq, Show)
newtype PValue = PValue {unPValue :: Double} deriving (Eq, Show)
newtype FETName = FETName {unFETName :: BS.ByteString} deriving(Eq, Show)
data FETMode = OneTail | TwoTail deriving (Eq,Show)

data FETResult = FETResult{groupOneA :: Int
                          ,groupOneB :: Int
                          ,groupTwoA :: Int
                          ,groupTwoB :: Int
                          ,pvalue    :: Double
                          ,fetName   :: BS.ByteString} deriving (Eq, Show)

fet :: FETName
    -> GroupOneA
    -> GroupOneB
    -> GroupTwoA
    -> GroupTwoB
    -> FETMode
    -> FETResult
fet (FETName fName) (GroupOneA goa) (GroupOneB gob) (GroupTwoA gta) (GroupTwoB gtb) mode
  | l > 0 && goa >= 0 && gob >= 0 && gta >= 0 && gtb >= 0 && k > 0 && gt > 0 =
    case mode of
        OneTail -> theResult{pvalue = cumD}
        TwoTail -> theResult{pvalue = twoTailValue}
  | otherwise = theResult
  where
    cumD = if l >= 1000
            then chiP / 2
            else cumulative d $ fromIntegral goa
    probX = probability d goa
    twoTailValue = if l >= 1000
                      then chiP
                      else sum . filter (<= probX) . fmap (probability d) $ [0..k]
    d = hypergeometric m l k
    chiP = complCumulative d' $ chiSquare goa gob gta gtb
    d' = chiSquared 1
    m = goa + gta
    l = m + gob + gtb
    k = goa + gob
    gt = gta + gtb
    theResult = FETResult{groupOneA = goa
                         ,groupOneB = gob
                         ,groupTwoA = gta
                         ,groupTwoB = gtb
                         ,pvalue = 1
                         ,fetName = fName
                         }

chiSquare :: Int
          -> Int
          -> Int
          -> Int
          -> Double
chiSquare a' b' c' d' = (x' * x' * (a + b + c + d)) /
                        ((a+b) * (c+d) * (b+d) * (a+c))
  where
    x' = (a * d) - (b * c)
    a = fromIntegral a'
    b = fromIntegral b'
    c = fromIntegral c'
    d = fromIntegral d'

max_stars_repo_name: chadlaing/feht

max_issues_repo_name: chadlaing/feht

max_forks_repo_name: chadlaing/feht


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Main where

import Lib
import Numeric.LinearAlgebra
import Data.Traversable

train_1_path = "/Users/klinkhae/OSU_Coursework/Rob537/HW1/hw1_data/train1.csv"
train_2_path = "/Users/klinkhae/OSU_Coursework/Rob537/HW1/hw1_data/train2.csv"
test_1_path = "/Users/klinkhae/OSU_Coursework/Rob537/HW1/hw1_data/test1.csv"
test_2_path = "/Users/klinkhae/OSU_Coursework/Rob537/HW1/hw1_data/test2.csv"
test_3_path = "/Users/klinkhae/OSU_Coursework/Rob537/HW1/hw1_data/test3.csv"

lrs = [0.005, 0.01, 0.05, 0.1]
hus = [4, 8, 12, 16]
eps = [1000, 5000, 10000, 100000]
sig = [True, False]

main :: IO ()
main = do
  let exps = generateExps lrs eps hus sig

  putStrLn "All Experiments with training data set 1"
  runExpsIO train_1_path [train_1_path, test_1_path, test_2_path, test_3_path] exps

  putStrLn "All Experiments with training data set 2"
  runExpsIO train_2_path [train_2_path, test_1_path, test_2_path, test_3_path] exps


text:
{-# LANGUAGE NumericUnderscores, OverloadedStrings #-}

module Main where

import Test.QuickCheck

import qualified Data.Vector.Unboxed as U
import qualified Data.Vector as V
-- import Data.Vector.Unboxed  ((!))
import Data.Time (fromGregorian, UTCTime)
import Data.Time.Format (parseTimeM, defaultTimeLocale)
import qualified Data.HashMap.Strict as Hm
import Data.Text (Text)
import qualified Data.Complex as Cx
import Data.Complex (Complex ( (:+) ) )

import Data.Approx

main :: IO ()
main = do

  -- print "Float equality test"
  quickCheck $ (0.0 :: Float) =~ (-1e-9)
  quickCheck $ (0.0 :: Double) =~ (-0.0)
  quickCheck $ (1.0e+7 :: Double) =~ 10_000_000.05

  quickCheck $ Just (10.0000000000000007 :: Double) =~ Just 10.0
  quickCheck $ Just (10 :: Double) /~ Just 1.0
  quickCheck $ Just (10 :: Double) /~ Nothing
  quickCheck $ (Nothing :: Maybe Double) /~ Just 1.0
  quickCheck $ (Nothing :: Maybe Double) =~ Nothing

  quickCheck $ ([1.2, 3.4, 5.6] :: [Double]) /~ [1.2, 3.4, 5.65]
  quickCheck $ ([2.3, 1,2, 1000] :: [Double]) =~ [2.3, 1,2, 1e+3]
  quickCheck $ ([1.0,1.0,1.0,1.0] :: [Double]) /~ [1.0,1.0,1.0,1.0,2.0]

  quickCheck $ ((1.2, 3.4) :: (Double, Double)) =~ (1.20000008, 3.399999999)
  quickCheck $ ((1.2, 3.5) :: (Double, Double)) /~ (1.20000008, 3.399999999)

  quickCheck $ ((1.2, 3.4, 2.5) :: (Double, Double, Double)) =~ (1.20000008, 3.399999999, 2.5000001)
  quickCheck $ ((1.2, 3.4, 2.5) :: (Double, Double, Double)) /~ (1.2, 3.399999999, 2.4)

  quickCheck $ ((1.2, 3.4, 2.5) :: (Float, Float, Float)) =~ (1.20000008, 3.399999999, 2.5000001)
  quickCheck $ ((1.2, 3.4, 2.5) :: (Float, Float, Float)) /~ (1.2, 3.399999999, 2.4)

  quickCheck $ ((1, 3, 2) :: (Int, Int, Int)) =~ (1, 3, 2)
  quickCheck $ ((1, 3, 2) :: (Int, Int, Int)) /~ (1, 3, 5)

  quickCheck $ ((2,5.35) :: (Int,Double)) =~ (2,5.35)

  quickCheck $ ((3.0 :+ 5.0)::(Cx.Complex Float)) =~ (3.0 :+ 5.0)

  quickCheck $ ((2,5.35,"happ",[1,2],["ret","we"],[6.78,3.5])
    :: (Int,Double,Text,[Int],[[Char]],[Double]))

    =~ (2,5.35,"happ",[1,2],["ret","we"],[6.78,3.5])

  quickCheck $ ((1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0) :: (Double,Double,Double,Double,Double,Double,Double,Double,Double,Double,Double)) =~ (1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0)

  quickCheck $ ((1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.1) :: (Double,Double,Double,Double,Double,Double,Double,Double,Double,Double,Double))   /~ (1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0,11.0)

  quickCheck $ ((1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0) :: (Double,Double,Double,Double,Double,Double,Double,Double,Double,Double)) =~ (1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0)

  quickCheck $ ((1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.0) :: (Double,Double,Double,Double,Double,Double,Double,Double,Double,Double))   /~ (1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10.1)

  quickCheck $ False /~ True
  quickCheck $ True =~ True
  quickCheck $ False =~ False

  quickCheck $ ("star"::Text) =~ ("star"::Text)
  quickCheck $ ("star"::Text) /~ ("start"::Text)
  quickCheck $ ("start"::Text) /~ ("star"::Text)

  quickCheck $ ("star"::String) =~ ("star"::String)
  quickCheck $ ("star"::String) /~ ("start"::String)
  quickCheck $ ("start"::String) /~ ("star"::String)

  quickCheck $ fromGregorian 2018 3 30 =~ fromGregorian 2018 3 30
  quickCheck $ fromGregorian 2018 3 31 /~ fromGregorian 2018 3 30

  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15" :: Maybe UTCTime) /~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 14:02:15"

  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15" :: Maybe UTCTime) =~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:01:50"

  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:01:50" :: Maybe UTCTime) =~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15"
  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-16 15:02:15" :: Maybe UTCTime) /~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15"

  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2019-01-15 15:02:15" :: Maybe UTCTime) /~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15"

  quickCheck $ (parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15" :: Maybe UTCTime) =~ parseTimeM True defaultTimeLocale "%Y-%m-%d %H:%M:%S" "2020-01-15 15:02:15"

  let rb = Hm.fromList [("Cash", 10.23), ("CurrentAdvances", 56.25), ("AccountPayables", 0.0)] :: Hm.HashMap String Double
  let tb = Hm.fromList [("Cash", 10.23), ("CurrentAdvances", 56.25)] :: Hm.HashMap String Double

  quickCheck $ rb /~ tb
  quickCheck $ tb /~ rb

  let rb = Hm.fromList [("Cash", 10.23), ("CurrentAdvances", 56.25)] :: Hm.HashMap String Double

  quickCheck $ rb =~ tb
  quickCheck $ tb =~ rb

  let rb = Hm.fromList [("Cash", 10.21), ("CurrentAdvances", 56.25)] :: Hm.HashMap String Double

  quickCheck $ rb /~ tb
  quickCheck $ tb /~ rb

  let rb = Hm.fromList [("Cash", Just 10.23), ("CurrentAdvances", Just 56.25), ("AccountPayables", Just 0.0)] :: Hm.HashMap String (Maybe Double)
  let tb = Hm.fromList [("Cash", Just 10.23), ("CurrentAdvances", Just 56.25)] :: Hm.HashMap String (Maybe Double)

  quickCheck $ rb /~ tb
  quickCheck $ tb /~ rb

  let rb = Hm.fromList [("Cash", Just 10.23), ("CurrentAdvances", Just 56.25)] :: Hm.HashMap String (Maybe Double)

  quickCheck $ rb =~ tb
  quickCheck $ tb =~ rb

  let rb = V.fromList [("Cash", Just 10.23), ("CurrentAdvances", Just 56.25), ("AccountPayables", Just 0.0)] :: V.Vector (String, Maybe Double)
  let tb = V.fromList [("Cash", Just 10.23), ("CurrentAdvances", Just 56.25)] :: V.Vector (String, Maybe Double)

  quickCheck $ rb /~ tb

  let rb = V.fromList [("Cash", Just 10.24), ("CurrentAdvances", Just 56.25)] :: V.Vector (String, Maybe Double)

  quickCheck $ rb /~ tb

  quickCheck $ U.fromList ([2.5, 1.5, 3.56] :: [Double]) =~ U.fromList [2.500000001, 1.5, 3.56]

  quickCheck $ (2,4) `inRange` (3::Int)
  quickCheck $ not $ (2,4) `inRange` (5::Int)

  quickCheck $ (2.5, 2.75) `inRange` (2.62::Double)
  quickCheck $ not $ (2.5, 2.75) `inRange` (2.8::Double)

  quickCheck $ (4,2) `safeInRange` (3::Int)
  quickCheck $ not $ (2,4) `safeInRange` (5::Int)

  quickCheck $ (2.5, 2.75) `safeInRange` (2.62::Double)
  quickCheck $ not $ (2.75, 2.5) `safeInRange` (2.8::Double)

  quickCheck $ inTol 0.25 6.25 (6.1::Double)
  quickCheck $ not $ inTol 0.125 6.25 (6.1::Double)

max_stars_repo_name: n-kishaloy/approx

max_issues_repo_name: n-kishaloy/approx

max_forks_repo_name: n-kishaloy/approx


text:
module SimpleR.Language.Syntax where

import qualified Data.Complex as C

type Complex = C.Complex Double

data IdentAnnot = IdentAnnot deriving (Ord, Eq, Show, Read)

data Ident = Ident
  { idName :: String
  , idPkg :: Maybe String
  , idAnnot :: Maybe IdentAnnot
  } deriving (Ord, Eq, Show, Read)

data TickAnnot = TickAnnot  deriving (Eq, Show, Read)

data Tick = Tick
  { tickAnnot :: TickAnnot
  }  deriving (Eq, Show, Read)

data Const =
    IntConst Int
  | DoubleConst Double
  | ComplexConst Complex
  | BoolConst Bool
  | StringConst String
  | NilConst
  | NaConst
   deriving (Eq, Show, Read)

data Param =
    Param Ident
  | Default Ident Expr
  | VarParam
   deriving (Eq, Show, Read)

data Arg =
    Arg Expr
  | Named Ident Expr
  | VarArg
   deriving (Eq, Show, Read)

data Expr =
    Var Ident
  | Const Const
  | Seq [Expr]
  | LamAbs [Param] Expr
  | LamApp Expr [Arg]
  | NativeLamApp Ident [Ident]
  | Return Expr
  | Assign Expr Expr

  | VecInd Expr [Arg]
  | VecIndAssign Expr [Arg] Expr
  | VecSub Expr [Arg]
  | VecSubAssign Expr [Arg] Expr
  | ListName Expr Expr
  | ListNameAssign Expr Expr Expr
  | ObjAttr Expr Expr
  | ObjAttrAssign Expr Expr Expr

  | If Expr Expr Expr
  | While Expr Expr
  | Break
  | Next

  | Error
  | Tix Tick Expr
   deriving (Eq, Show, Read)

data Program = Program [Expr]  deriving (Eq, Show, Read)

max_stars_repo_name: AntonXue/simple-r

max_issues_repo_name: AntonXue/simple-r

max_forks_repo_name: AntonXue/simple-r


text:
{-# LANGUAGE ScopedTypeVariables #-}
module Language.Haskell.Repl
    ( Repl
    -- * Construction
    , newRepl
    , repl'
    , defaultFlags
    , defaultImports
    , defaultLineLength
    , defaultPatience
    , defaultBuildExpr
    , defaultProcessOutput
    -- * Stopping
    , stopRepl
    -- * Interaction
    , Input(..)
    , ReplOutput(..)
    , Output(..)
    , prompt
    , prompt_
    , send
    , result
    , prettyOutput
    , parseInput
    ) where

import Control.Concurrent
import Control.Applicative ((<$>))
import Control.Exception (catch, SomeException(..), ErrorCall(..), fromException, Exception(..), evaluate)
import Control.Monad
import Control.Arrow
import Data.Dynamic
import Data.IORef
import Data.Char (isAscii)
import Data.Maybe
import Text.Parsec hiding (newline)
import Text.Parsec.String
import qualified Language.Haskell.Exts as H
import GHC
import GHC.Paths
import GhcMonad
import Outputable (showSDocForUser, Outputable, ppr, neverQualify)
import HscTypes
import OccName
import System.IO.Unsafe

data Repl a = Repl
    { inputChan         :: Chan Input
    , outputChan        :: Chan (ReplOutput a)
    , interpreter       :: ThreadId
    , config            :: Config a
    }

data Config a = Config
    { imports           :: [String]
    , compilerFlags     :: [String]
    , buildExpr         :: String  -> String
    , process           :: Dynamic -> IO a
    , patience          :: Double
    , lineLength        :: Int
    }

data Input
    = Type     String
    | Kind     String
    | Info     String
    | Decl     String
    | Stmt     String
    | Expr     String
    | Undefine String
    | Clear
  deriving Show

data ReplOutput a
    = ReplError String
    | GhcError  String
    | Output    [String]
    | Result    a
  deriving Show

data Output
    = OK        [String]
    | Exception [String] String
    | Errors    [String]
    | Partial   [String]
    | Timeout   [String]
  deriving Show

prefix :: String -> Parser ()
prefix (x:xs) = do
    _ <- string [':',x]
    forM_ xs (optional . char)
    spaces
prefix [] = fail "empty prefix"

simpl :: String -> (String -> a) -> Parser a
simpl pfix f = do
    prefix pfix
    f <$> getInput

parseType, parseKind, parseInfo, parseDecl, parseStmt, parseExpr, parseUndefine, parseClear, parseInput :: Parser Input
parseType     = simpl "type"  Type
parseKind     = simpl "kind"  Kind
parseInfo     = simpl "info"  Info
parseUndefine = simpl "undef" Undefine
parseClear    = simpl "clear" (const Clear)
parseDecl     = do
    -- from InteractiveUI.hs
    p <- single ["class ","type ","data ","newtype ","instance ","deriving ","foreign ","default(","default "]
    r <- getInput
    return (Decl (p ++ r))
  where single = foldr1 (<|>) . map (try . string)

parseStmt = do
    -- Problem: a Stmt is automatically ran if it is :: IO a
    -- So we have to make sure it is a let binding.
    -- BUT haskell-src-exts can't handle Unicode in let bindings, so valid
    -- let bindings like "let あ = 0" get obliterated. Therefore, short of
    -- actually parsing the binding ourselves, we replaceanything not ASCII
    -- with an ASCII character ('x').
    stmt <- getInput
    case H.parseStmt stmt of
        H.ParseOk H.LetStmt{}
          -> return (Stmt stmt)
        _ -> case H.parseStmt (mangle stmt) of
            H.ParseOk H.LetStmt{}
              -> return (Stmt stmt)
            _ -> fail "Not a let binding."
  where
    mangle = map $ \c -> if isAscii c then c else 'x'

parseExpr = Expr <$> getInput

-- | Used by 'prompt'
parseInput   = foldr1 (\l r -> Text.Parsec.try l <|> r)
    [ parseClear
    , parseUndefine
    , parseType
    , parseKind
    , parseInfo
    , parseStmt
    , parseDecl
    , parseExpr ]

unsafeCatch :: Exception e => a -> (e -> a) -> a
unsafeCatch a f = unsafePerformIO (catch (evaluate a) (return . f))

cripple :: a -> a -> a
cripple x y = unsafeCatch x (\SomeException{} -> y)

-- | Used by 'prompt'.
prettyOutput :: Repl a -> Output -> [String]
prettyOutput _ (OK s)          = s
prettyOutput _ (Partial s)     = s
prettyOutput _ (Errors errs)   = errs
prettyOutput r (Exception s e) = map
    (take (lineLength (config r)))
    (overLast (++ ("*** Exception: " ++ cripple e "*** Exception: that's enough exceptions for you.")) s)
prettyOutput _ (Timeout [])    = ["*** Timed out"]
prettyOutput _ (Timeout s)     = overLast (++ "*** Timed out") s

-- | Send input.
send :: Repl a -> Input -> IO ()
send = writeChan . inputChan

-- | Naiively get the next set of results. This /does not/ take into account
-- 'patiences', 'patienceForErrors', or 'lineLength'. However, due
-- to laziness, this may not matter.
result :: Repl a -> IO (ReplOutput a)
result = readChan . outputChan

{-# INLINE index #-}
index :: Int -> [a] -> Maybe a
i `index` ys
    | i >= 0    = go 0 ys
    | otherwise = Nothing
  where
    go _ [] = Nothing
    go j (x:xs)
        | j == i    = Just x
        | otherwise = go (j+1) xs

{-# INLINE overHead #-}
overHead :: (a -> a) -> [a] -> [a]
overHead f xs' = case xs' of
    x:xs -> f x : xs
    _    -> []

{-# INLINE overLast #-}
overLast :: (a -> a) -> [a] -> [a]
overLast f = go
  where
    go []     = []
    go [x]    = [f x]
    go (x:xs) = x : go xs

{-# INLINE lengthAt #-}
lengthAt :: Int -> [[a]] -> Int
lengthAt i = maybe 0 length . index i

-- | Same as 'prompt_', except it parses the input, and pretty prints the results.
prompt
    :: Repl [String]
    -> String
    -> IO [String]
prompt repl x = prettyOutput repl <$> prompt_ repl (case runParser parseInput () "" x of
    Right a -> a
    -- Should be impossible to reach. parseExpr gobbles up everything.
    _       -> error "Cannot parse input!")

-- | Enter commands and an expression to a 'Repl', and immediately consume results.
-- However, truncate input to the given length, and stop the computation after the
-- given amount of time in seconds.
prompt_
    :: Repl [String]
    -> Input
    -> IO Output
prompt_ repl x = do
    send repl x
    results' <- result repl

    unlessRedundant results' $ \ results -> do
        -- outputs is used iff an exception is raised by the compiled input.
        -- This was a DList, but I didn't find any real advantage of it over
        -- [String] -- snoc was cheap but toList very expensive.
        outputs :: IORef [String]   <- newIORef []
        threads :: IORef [ThreadId] <- newIORef []
        final   :: MVar Output      <- newEmptyMVar
        let push    c   = do
                output <- readIORef outputs
                if lengthAt (length output - 1) output > lineLength (config repl)
                  then putMVar final (Partial (unreverse output))
                  else writeIORef outputs (overHead (c:) output)
            newline     = modifyIORef outputs ([]:)
            readOutput  = unreverse <$> readIORef outputs
            fork f      = do
                thread <- forkIO $ f `catch` \e@SomeException{} -> do
                    output <- readOutput
                    putMVar final (Exception output (show e))
                modifyIORef threads (thread:)

        -- Time out
        -- This can return only Timeout <what was consumed so far>
        fork $ do
            threadDelay (floor (patience (config repl) * 1000000))
            output <- readOutput
            putMVar final (Timeout output)

        -- Read characters off of results, and "push" them to outputs.
        fork $ do
            forM_ results $ \l -> do
                newline
                forM_ l push
            putMVar final . OK =<< readOutput

        output <- takeMVar final
        mapM_ killThread =<< readIORef threads
        return output
  where
    unreverse = reverse . map reverse
    trim      = take (lineLength (config repl))

    -- | Don't bother with things other than an actual result from an expression -- they will be loaded "instantly"
    unlessRedundant (ReplError s) _ = return . Errors . map trim . lines $ s
    unlessRedundant (GhcError  s) _ = return . Errors . map trim . lines $ s
    unlessRedundant (Output s) _    = return . OK . map trim $ s
    unlessRedundant (Result s) f    = f s

stopRepl :: Repl a -> IO ()
stopRepl = killThread . interpreter

-- | Default GHCi-like 'Repl'
newRepl :: IO (Repl [String])
newRepl = do
    inp <- newChan
    out <- newChan
    repl' inp out defaultConfig

defaultConfig :: Config [String]
defaultConfig = Config
    defaultImports
    defaultFlags
    defaultBuildExpr
    defaultProcessOutput
    defaultPatience
    defaultLineLength

defaultImports :: [String]
defaultImports
  = ["import Prelude hiding ((.), id, catch)"
    ,"import GHC.TypeLits"
    ,"import qualified Data.Map as M"
    ,"import qualified Data.Foldable as F"
    ,"import qualified Data.Traversable as T"
    ,"import qualified Control.Exception (catch)"
    ,"import Control.Monad.Reader"
    ,"import Control.Monad.State"
    ,"import Control.Monad.Writer"
    ,"import Control.Monad.RWS"
    ,"import Control.Monad.Identity"
    ,"import Control.Monad.ST"
    ,"import Control.Comonad"
    ,"import Control.Category"
    ,"import Control.Monad"
    ,"import Control.Monad.Fix"
    ,"import Control.Applicative"
    ,"import Control.Lens"
    ,"import Control.Arrow"
    ,"import Data.Function hiding ((.), id)"
    ,"import Data.Either"
    ,"import Data.Int"
    ,"import Data.Word"
    ,"import Data.List"
    ,"import Data.List.Split"
    ,"import Data.Maybe"
    ,"import Data.Bits"
    ,"import Data.Array"
    ,"import Data.Ix"
    ,"import Data.Functor"
    ,"import Data.Typeable"
    ,"import Data.Monoid"
    ,"import Data.Ratio"
    ,"import Data.Complex"
    ,"import Data.Char"
    ,"import Data.Bits.Lens"
    ,"import Data.List.Lens"
    ,"import Data.List.Split.Lens"
    ]

defaultFlags :: [String]
defaultFlags = map ("-X"++)
    ["DataKinds"
    ,"PolyKinds"
    ,"KindSignatures"
    ,"TypeOperators"
    ,"DeriveFunctor"
    ,"DeriveTraversable"
    ,"DeriveFoldable"
    ,"DeriveDataTypeable"
    ,"DeriveGeneric"
    ,"OverloadedStrings"
    ,"ImplicitParams"
    ,"BangPatterns"
    ,"PatternGuards"
    ,"MultiWayIf"
    ,"LambdaCase"
    ,"FlexibleInstances"
    ,"FlexibleContexts"
    ,"FunctionalDependencies"
    ,"StandaloneDeriving"
    ,"MultiParamTypeClasses"
    ,"UnicodeSyntax"
    ,"RankNTypes"
    ,"ExistentialQuantification"
    ,"GADTs"
    ,"TypeFamilies"
    ,"Safe"
    ] ++
    [ "-dcore-lint" ]

-- | defaultLineLength = 512
defaultLineLength :: Int
defaultLineLength = 512

-- | defaultPatience = 5
defaultPatience :: Double
defaultPatience = 5

defaultBuildExpr :: String -> String
defaultBuildExpr x = "show (" ++ x ++ ")"

defaultProcessOutput :: Dynamic -> IO [String]
defaultProcessOutput d = return (lines (fromDyn d ""))

-- | 'Repl' smart constructor.
repl'
    :: Chan Input           -- ^ Input channel
    -> Chan (ReplOutput a)  -- ^ Output channel
    -> Config a
    -> IO (Repl a)
repl' inp out cfg = do
    interp <- forkIO $
        runGhc (Just libdir) $ do
            initialDynFlags <- getProgramDynFlags
            (dflags',_,_)   <- parseDynamicFlags initialDynFlags (map (mkGeneralLocated "flag") (compilerFlags cfg))
            _pkgs           <- setSessionDynFlags dflags'
            dflags          <- getSessionDynFlags

            let sdoc :: Outputable a => a -> String
                sdoc = showSDocForUser dflags neverQualify . ppr

                formatType
                    = splitForAllTys
                  >>> snd
                  >>> sdoc
                  >>> lines
                  >>> Output

            forever $ do
                import_ (imports cfg)
                i'   <- liftIO (readChan inp)
                liftIO . writeChan out =<< case i' of
                    Clear      -> do
                        setTargets []
                        void (load LoadAllTargets)
                        return (Output ["Cleared memory."])

                    Undefine s' -> fmap Output $
                        forM (words s') $ \s -> do
                            let eqs :: NamedThing a => a -> Bool
                                eqs n = occNameString (getOccName n) == s
                            session <- getSession
                            setSession session
                                { hsc_IC = (hsc_IC session)
                                    { ic_tythings = filter (not . eqs) (ic_tythings (hsc_IC session)) }
                                }
                            return "OK."

                    Type s -> errors $ formatType <$> exprType s
                    Kind s -> errors $ formatType . snd <$> typeKind True s
                    Decl s -> errors $ do _names <- runDecls s; return (Output ["OK."]) --  ["OK."]
                    Stmt s -> errors $ do void (runStmt s SingleStep); return (Output ["OK."])
                    Expr s -> errors $ do
                        compiled <- dynCompileExpr (buildExpr cfg s)
                        built    <- liftIO (process cfg compiled)
                        return (Result built)
                    Info s -> errors $ do
                        names <- parseName s
                        infos <- concatMap (\(t,f,cs) -> sdoc t : sdoc f : map sdoc cs)
                               . catMaybes
                             <$> mapM getInfo names
                        return $ Output infos
    return Repl
        { inputChan      = inp
        , outputChan     = out
        , interpreter    = interp
        , config         = cfg
        }
  where
    errors x = x `gcatch` \ e@SomeException{} ->
        return $! case fromException e :: Maybe ErrorCall of
            Just _ -> ReplError (show e)
            _      -> GhcError  (show e)

    import_ = mapM (fmap IIDecl . parseImportDecl) >=> setContext


text:
{-# LANGUAGE FlexibleContexts #-}

module Evaluator.Numerical where

import LispTypes
import Environment
import Evaluator.Operators

import Data.Complex
import Data.Ratio
import Data.Foldable
import Data.Fixed
import Numeric
import Control.Monad.Except

numericalPrimitives :: [(String, [LispVal] -> ThrowsError LispVal)]
numericalPrimitives =
    -- Binary Numerical operations
    [("+", numAdd),
    ("-", numSub),
    ("*", numMul),
    ("/", numDivide),
    ("modulo", numMod),
    ("quotient", numericBinop quot),
    ("remainder", numericBinop rem),
    ("abs", numAbs),
    ("ceiling", numCeil),
    ("floor", numFloor),
    ("round", numRound),
    ("truncate", numTruncate),
    -- Conversion
    ("number->string", numToString),
    -- Numerical Boolean operators
    ("=", numBoolBinop (==)),
    ("<", numBoolBinop (<)),
    (">", numBoolBinop (>)),
    ("/=", numBoolBinop (/=)),
    (">=", numBoolBinop (>=)),
    ("<=", numBoolBinop (<=)),
    -- Scientific functions
    ("sqrt", sciSqrt),
    ("acos", sciAcos),
    ("asin", sciAsin),
    ("atan", sciAtan),
    ("cos", sciCos),
    ("sin", sciSin),
    ("tan", sciTan),
--    ("acosh", sciAcosh),
--    ("asinh", sciAsinh),
--    ("atanh", sciAtanh),
--    ("cosh", sciCosh),
--    ("sinh", sciSinh),
--    ("tanh", sciTanh),
    ("exp", sciExp),
    ("expt", sciExpt),
    ("log", sciLog),
    -- Complex numbers operations
--    ("angle", cAngle),
    ("real-part", cRealPart),
    ("imag-part", cImagPart),
    ("magnitude", cMagnitude),
    ("make-polar", cMakePolar),
--    ("make-rectangular", cMakeRectangular),
    -- Type testing functions
    ("number?", unaryOp numberp),
    ("integer?", unaryOp integerp),
    ("float?", unaryOp floatp),
    ("ratio?", unaryOp ratiop),
    ("complex?", unaryOp complexp)]

-- |Type testing functions
numberp, integerp, floatp, ratiop, complexp :: LispVal -> LispVal
integerp (Number _)     = Bool True
integerp _              = Bool False
numberp (Number _)      = Bool True
numberp (Float _)       = Bool True
numberp (Ratio _)       = Bool True
numberp (Complex _)     = Bool True
numberp _               = Bool False
floatp (Float _)        = Bool True
floatp _                = Bool False
ratiop (Ratio _)        = Bool True
ratiop _                = Bool False
complexp (Complex _)    = Bool True
complexp _              = Bool False

-- |Absolute value
numAbs :: [LispVal] -> ThrowsError LispVal
numAbs [Number x]   = return $ Number $ abs x
numAbs [Float x]    = return $ Float $ abs x
numAbs [Complex x]  = return $ Complex $ abs x -- Calculates the magnitude
numAbs [Ratio x]    = return $ Ratio $ abs x
numAbs [x] = throwError $ TypeMismatch "number" x
numAbs l = throwError $ NumArgs 1 l

-- |Ceiling, floor, round and truncate
numCeil, numFloor, numRound, numTruncate :: [LispVal] -> ThrowsError LispVal
numCeil [Number x]   = return $ Number $ ceiling $ fromInteger x
numCeil [Ratio x]    = return $ Number $ ceiling $ fromRational x
numCeil [Float x]    = return $ Number $ ceiling x
numCeil [Complex x]  =
    if imagPart x == 0 then return $ Number $ ceiling $ realPart x
    else throwError $ TypeMismatch "integer or float" $ Complex x
numCeil [x] = throwError $ TypeMismatch "integer or float" x
numCeil l = throwError $ NumArgs 1 l

numFloor [Number x]   = return $ Number $ floor $ fromInteger x
numFloor [Ratio x]    = return $ Number $ floor $ fromRational x
numFloor [Float x]    = return $ Number $ floor x
numFloor [Complex x]  =
    if imagPart x == 0 then return $ Number $ floor $ realPart x
    else throwError $ TypeMismatch "integer or float" $ Complex x
numFloor [x] = throwError $ TypeMismatch "integer or float" x
numFloor l = throwError $ NumArgs 1 l

numRound [Number x]   = return $ Number $ round $ fromInteger x
numRound [Ratio x]    = return $ Number $ round $ fromRational x
numRound [Float x]    = return $ Number $ round x
numRound [Complex x]  =
    if imagPart x == 0 then return $ Number $ round $ realPart x
    else throwError $ TypeMismatch "integer or float" $ Complex x
numRound [x] = throwError $ TypeMismatch "integer or float" x
numRound l = throwError $ NumArgs 1 l

numTruncate [Number x]   = return $ Number $ truncate $ fromInteger x
numTruncate [Ratio x]    = return $ Number $ truncate $ fromRational x
numTruncate [Float x]    = return $ Number $ truncate x
numTruncate [Complex x]  =
    if imagPart x == 0 then return $ Number $ truncate $ realPart x
    else throwError $ TypeMismatch "integer or float" $ Complex x
numTruncate [x] = throwError $ TypeMismatch "integer or float" x
numTruncate l = throwError $ NumArgs 1 l

-- |foldl1M is like foldlM but has no base case
foldl1M :: Monad m => (a -> a -> m a) -> [a] -> m a
foldl1M f (x : xs) = foldlM f x xs
foldl1M _ _ = error "unexpected error in foldl1M"

-- | Sum numbers
numAdd :: [LispVal] -> ThrowsError LispVal
numAdd []   = return $ Number 0
numAdd l    = foldl1M (\ x y -> numCast [x, y] >>= go) l
    where
        go (List [Number x, Number y])      = return $ Number $ x + y
        go (List [Float x, Float y])        = return $ Float $ x + y
        go (List [Ratio x, Ratio y])        = return $ Ratio $ x + y
        go (List [Complex x, Complex y])    = return $ Complex $ x + y
        go _ = throwError $ Default "unexpected error in (+)"

-- | Subtract numbers
numSub :: [LispVal] -> ThrowsError LispVal
numSub []           = throwError $ NumArgs 1 []
numSub [Number x]   = return $ Number $ -1 * x
numSub [Float x]    = return $ Float $ -1 * x
numSub [Ratio x]    = return $ Ratio $ -1 * x
numSub [Complex x]  = return $ Complex $ -1 * x
numSub l = foldl1M (\ x y -> numCast [x, y] >>= go) l
    where
        go (List [Number x, Number y])      = return $ Number $ x - y
        go (List [Float x, Float y])        = return $ Float $ x - y
        go (List [Ratio x, Ratio y])        = return $ Ratio $ x - y
        go (List [Complex x, Complex y])    = return $ Complex $ x - y
        go _ = throwError $ Default "unexpected error in (-)"

-- | Multiply numbers
numMul :: [LispVal] -> ThrowsError LispVal
numMul []   = return $ Number 1
numMul l    = foldl1M (\ x y -> numCast [x, y] >>= go) l
    where
        go (List [Number x, Number y])      = return $ Number $ x * y
        go (List [Float x, Float y])        = return $ Float $ x * y
        go (List [Ratio x, Ratio y])        = return $ Ratio $ x * y
        go (List [Complex x, Complex y])    = return $ Complex $ x * y
        go _ = throwError $ Default "unexpected error in (*)"

-- |Divide two numbers
numDivide :: [LispVal] -> ThrowsError LispVal
numDivide [] = throwError $ NumArgs 1 []
numDivide [Number 0]   = throwError DivideByZero
numDivide [Ratio 0]     = throwError DivideByZero
numDivide [Number x]   = return $ Ratio $ 1 / fromInteger x
numDivide [Float x]    = return $ Float $ 1.0 / x
numDivide [Ratio x]    = return $ Ratio $ 1 / x
numDivide [Complex x]  = return $ Complex $ 1 / x
numDivide l = foldl1M (\ x y -> numCast [x, y] >>= go) l
    where
        go (List [Number x, Number y])
            | y == 0 = throwError DivideByZero
            | mod x y == 0 = return $ Number $ div x y -- Integer division
            | otherwise = return $ Ratio $ fromInteger x / fromInteger y
        go (List [Float x, Float y])
            | y == 0 = throwError DivideByZero
            | otherwise = return $ Float $ x / y
        go (List [Ratio x, Ratio y])
            | y == 0 = throwError DivideByZero
            | otherwise = return $ Ratio $ x / y
        go (List [Complex x, Complex y])
            | y == 0 = throwError DivideByZero
            | otherwise = return $ Complex $ x / y
        go _ = throwError $ Default "unexpected error in (/)"

-- |Numerical modulus
numMod :: [LispVal] -> ThrowsError LispVal
numMod [] = return $ Number 1
numMod l = foldl1M (\ x y -> numCast [x, y] >>= go) l
    where
        go (List [Number a, Number b]) = return $ Number $ mod' a b
        go (List [Float a, Float b]) = return $ Float $ mod' a b
        go (List [Ratio a, Ratio b]) = return $ Ratio $ mod' a b
        -- #TODO implement modulus for complex numbers
        go (List [Complex a, Complex b]) = throwError $ Default "modulus is not yet implemented for complex numbers"
        go _ = throwError $ Default "unexpected error in (modulus)"

-- | Boolean operator
numBoolBinop :: (LispVal -> LispVal -> Bool) -> [LispVal] -> ThrowsError LispVal
numBoolBinop op [] = throwError $ Default "need at least two arguments"
numBoolBinop op [x] = throwError $ Default "need at least two arguments"
numBoolBinop op [x, y] = numCast [x, y] >>= go op
    where
        go op (List [x@(Number _),  y@(Number _)])      = return $ Bool $ x `op` y
        go op (List [x@(Float _),   y@(Float _)])       = return $ Bool $ x `op` y
        go op (List [x@(Ratio _),   y@(Ratio _)])       = return $ Bool $ x `op` y
        go op (List [x@(Complex _), y@(Complex _)])     = return $ Bool $ x `op` y
        go op _ = throwError $ Default "unexpected error in boolean operation"

-- |Accept two numbers and cast one as the appropriate type
numCast :: [LispVal] -> ThrowsError LispVal
-- Same type, just return the two numbers
numCast [a@(Number _),  b@(Number _)]   = return $ List [a,b]
numCast [a@(Float _),   b@(Float _)]    = return $ List [a,b]
numCast [a@(Ratio _),   b@(Ratio _)]    = return $ List [a,b]
numCast [a@(Complex _), b@(Complex _)]  = return $ List [a,b]
-- First number is an integer
numCast [Number a,  b@(Float _)]      = return $ List [Float $ fromInteger a, b]
numCast [Number a,  b@(Ratio _)]      = return $ List [Ratio $ fromInteger a, b]
numCast [Number a,  b@(Complex _)]    = return $ List [Complex $ fromInteger a, b]
-- First number is a float
numCast [a@(Float _),  Number b]      = return $ List [a, Float $ fromInteger b]
numCast [a@(Float _),  Ratio b]       = return $ List [a, Float $ fromRational b]
numCast [Float a,  b@(Complex _)]     = return $ List [Complex $ a :+ 0, b]
-- First number is a rational
numCast [a@(Ratio _),  Number b]      = return $ List [a, Ratio $ fromInteger b]
numCast [Ratio a,  b@(Float _)]       = return $ List [Float $ fromRational a, b]
numCast [Ratio a,  b@(Complex _)] =
    return $ List [Complex $ fromInteger (numerator a) / fromInteger (denominator a), b]
-- First number is a complex
numCast [a@(Complex _),  Number b] = return $ List [a, Complex $ fromInteger b]
numCast [a@(Complex _),  Float b]     = return $ List [a, Complex $ b :+ 0]
numCast [a@(Complex _),  Ratio b] =
    return $ List [a, Complex $ fromInteger (numerator b) / fromInteger (denominator b)]
-- Error cases
numCast [a, b] = case a of
        Number _ -> throwError $ TypeMismatch "number" b
        Float _ -> throwError $ TypeMismatch "number" b
        Ratio _ -> throwError $ TypeMismatch "number" b
        Complex _ -> throwError $ TypeMismatch "number" b
        _ -> throwError $ TypeMismatch "number" a
numCast _ = throwError $ Default "unknown error in numCast"

-- |Take a primitive Haskell Integer function and wrap it
-- with code to unpack an argument list, apply the function to it
-- and return a numeric value
numericBinop :: (Integer -> Integer -> Integer) -> [LispVal] -> ThrowsError LispVal
-- Throw an error if there's only one argument
numericBinop op val@[_] = throwError $ NumArgs 2 val
-- Fold the operator leftway if there are enough args
numericBinop op params = mapM unpackNum params >>= return . Number . foldl1 op

-- |Convert a Number to a String
numToString :: [LispVal] -> ThrowsError LispVal
numToString [Number x] = return $ String $ show x
numToString [Float x] = return $ String $ show x
numToString [Ratio x] = return $ String $ show x
numToString [Complex x] = return $ String $ show x
numToString [x] = throwError $ TypeMismatch "number" x
numToString badArgList = throwError $ NumArgs 2 badArgList

-- | Trigonometric functions
sciCos, sciSin, sciTan, sciAcos, sciAsin, sciAtan :: [LispVal] -> ThrowsError LispVal
-- | Cosine of a number
sciCos [Number x] = return $ Float $ cos $ fromInteger x
sciCos [Float x] = return $ Float $ cos x
sciCos [Ratio x] = return $ Float $ cos $ fromRational x
sciCos [Complex x] = return $ Complex $ cos x
sciCos [notnum] = throwError $ TypeMismatch "number" notnum
sciCos badArgList = throwError $ NumArgs 1 badArgList
-- | Sine of a number
sciSin [Number x] = return $ Float $ sin $ fromInteger x
sciSin [Float x] = return $ Float $ sin x
sciSin [Ratio x] = return $ Float $ sin $ fromRational x
sciSin [Complex x] = return $ Complex $ sin x
sciSin [notnum] = throwError $ TypeMismatch "number" notnum
sciSin badArgList = throwError $ NumArgs 1 badArgList
-- | Tangent of a number
sciTan [Number x] = return $ Float $ tan $ fromInteger x
sciTan [Float x] = return $ Float $ tan x
sciTan [Ratio x] = return $ Float $ tan $ fromRational x
sciTan [Complex x] = return $ Complex $ tan x
sciTan [notnum] = throwError $ TypeMismatch "number" notnum
sciTan badArgList = throwError $ NumArgs 1 badArgList
-- | Arccosine of a number
sciAcos [Number x] = return $ Float $ acos $ fromInteger x
sciAcos [Float x] = return $ Float $ acos x
sciAcos [Ratio x] = return $ Float $ acos $ fromRational x
sciAcos [Complex x] = return $ Complex $ acos x
sciAcos [notnum] = throwError $ TypeMismatch "number" notnum
sciAcos badArgList = throwError $ NumArgs 1 badArgList
-- | Sine of a number
sciAsin [Number x] = return $ Float $ asin $ fromInteger x
sciAsin [Float x] = return $ Float $ asin x
sciAsin [Ratio x] = return $ Float $ asin $ fromRational x
sciAsin [Complex x] = return $ Complex $ asin x
sciAsin [notnum] = throwError $ TypeMismatch "number" notnum
sciAsin badArgList = throwError $ NumArgs 1 badArgList
-- | Tangent of a number
sciAtan [Number x] = return $ Float $ atan $ fromInteger x
sciAtan [Float x] = return $ Float $ atan x
sciAtan [Ratio x] = return $ Float $ atan $ fromRational x
sciAtan [Complex x] = return $ Complex $ atan x
sciAtan [notnum] = throwError $ TypeMismatch "number" notnum
sciAtan badArgList = throwError $ NumArgs 1 badArgList

-- #TODO Implement hyperbolic functions
-- Ask teacher
-- | Hyperbolic functions
-- sciAcosh, sciAsinh, sciAtanh, sciCosh, sciSinh, sciTanh :: [LispVal] -> ThrowsError LispVal

-- Misc. Scientific primitives
sciSqrt, sciExp, sciExpt, sciLog :: [LispVal] -> ThrowsError LispVal
-- | Square root of a number
sciSqrt [Number x] = return $ Float $ sqrt $ fromInteger x
sciSqrt [Float x] = return $ Float $ sqrt x
sciSqrt [Ratio x] = return $ Float $ sqrt $ fromRational x
sciSqrt [Complex x] = return $ Complex $ sqrt x
sciSqrt [notnum] = throwError $ TypeMismatch "number" notnum
sciSqrt badArgList = throwError $ NumArgs 1 badArgList
-- | Return e to the power of x
sciExp [Number x] = return $ Float $ exp $ fromInteger x
sciExp [Float x] = return $ Float $ exp x
sciExp [Ratio x] = return $ Float $ exp $ fromRational x
sciExp [Complex x] = return $ Complex $ exp x
sciExp [notnum] = throwError $ TypeMismatch "number" notnum
sciExp badArgList = throwError $ NumArgs 1 badArgList
-- | Return x to the power of y
sciExpt [x, y] = numCast [x, y] >>= go
    where
        go (List [Number x, Number y])      = return $ Number $ round $ (fromInteger x) ** (fromInteger y)
        go (List [Float x, Float y])        = return $ Float $ x ** y
        go (List [Ratio x, Ratio y])        = return $ Float $ (fromRational x) ** (fromRational y)
        go (List [Complex x, Complex y])    = return $ Complex $ x ** y
        go _ = throwError $ Default "unexpected error in (-)"
sciExpt badArgList = throwError $ NumArgs 2 badArgList
-- | Return the natural logarithm of x
sciLog [Number x] = return $ Float $ log $ fromInteger x
sciLog [Float x] = return $ Float $ log x
sciLog [Ratio x] = return $ Float $ log $ fromRational x
sciLog [Complex x] = return $ Complex $ log x
sciLog [notnum] = throwError $ TypeMismatch "number" notnum
sciLog badArgList = throwError $ NumArgs 1 badArgList

-- #TODO implement phase (angle)
-- Ask teacher how to convert phase formats from haskell to guile scheme
-- | Complex number functions
cRealPart, cImagPart, cMakePolar, cMagnitude  :: [LispVal] -> ThrowsError LispVal
-- | Real part of a complex number
cRealPart [Number x] = return $ Number $ fromInteger x
cRealPart [Float x] = return $ Float x
cRealPart [Ratio x] = return $ Float $ fromRational x
cRealPart [Complex x] = return $ Float $ realPart x
cRealPart [notnum] = throwError $ TypeMismatch "number" notnum
cRealPart badArgList = throwError $ NumArgs 1 badArgList
-- | Imaginary part of a complex number
cImagPart [Number x] = return $ Number 0
cImagPart [Float x] = return $ Number 0
cImagPart [Ratio x] = return $ Number 0
cImagPart [Complex x] = return $ Float $ imagPart x
cImagPart [notnum] = throwError $ TypeMismatch "number" notnum
cImagPart badArgList = throwError $ NumArgs 1 badArgList
-- | Form a complex number from polar components of magnitude and phase.
cMakePolar [mag, p] = numCast [mag, p] >>= go
        where
            go (List [Number mag, Number p])
                | mag == 0 = return $ Number 0
                | p == 0 = return $ Number mag
                | otherwise = return $ Complex $ mkPolar (fromInteger mag) (fromInteger p)
            go (List [Float mag, Float p])
                | mag == 0 = return $ Number 0
                | p == 0 = return $ Float mag
                | otherwise = return $ Complex $ mkPolar mag p
            go (List [Ratio mag, Ratio p])
                | mag == 0 = return $ Number 0
                | p == 0 = return $ Float (fromRational mag)
                | otherwise = return $ Complex $ mkPolar (fromRational mag) (fromRational p)
            go val@(List [Complex mag, Complex p]) = throwError $ TypeMismatch "real" val
            go _ = throwError $ Default "unexpected error in make-polar"
cMakePolar badArgList = throwError $ NumArgs 2 badArgList
-- | Return the magnitude (length) of a complex number
cMagnitude [Number x] = return $ Number $ fromInteger x
cMagnitude [Float x] = return $ Float x
cMagnitude [Ratio x] = return $ Float $ fromRational x
cMagnitude [Complex x] = return $ Float $ magnitude x
cMagnitude [notnum] = throwError $ TypeMismatch "number" notnum
cMagnitude badArgList = throwError $ NumArgs 1 badArgList

max_stars_repo_name: zfnmxt/yasih

max_issues_repo_name: zfnmxt/yasih

max_forks_repo_name: zfnmxt/yasih


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveFoldable #-}
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}

module PLC1 (PLC1(..)) where

import Control.Applicative
import Control.Comonad
import Control.Comonad.Store
import Control.Exception (assert)
-- import Control.Monad
-- import Data.Distributive
import Data.Proxy
import Data.Unfoldable
import Data.Validity
import qualified Data.Vector as V
import Data.VectorSpace
import GHC.TypeLits
-- import Linear.Metric
-- import Linear.Vector
import qualified Numeric.LinearAlgebra as L
import qualified Test.QuickCheck as QC
import qualified Test.QuickCheck.Function as QC

import ComonadicVector
import Function
import Integration

-- https://hackage.haskell.org/package/polynomial-0.7.3

newtype PLC1 (n' :: Nat) c s = PLC1
  { getPLC1 :: ComonadicVector V.Vector s
  } deriving (Eq, Ord, Read, Show, Foldable, Functor, Traversable)

instance (KnownNat n', Validity s) => Validity (PLC1 n' c s) where
    validate (PLC1 xs) =
        mconcat
        [ n > 0 <?@> "length is positive"
        , length xs == n <?@> "vector has correct length "
        , xs <?!> "vector"
        ]
        where n = fromInteger (natVal (Proxy :: Proxy n'))
    isValid = isValidByValidating

instance KnownNat n' => Unfoldable (PLC1 n' c) where
    -- unfold xs = PLC1 <$> unfold xs
    -- unfold xs = (PLC1 . ComonadicVector 0 . V.fromListN n) <$> mklist n
    --     where n = fromInteger (natVal (Proxy :: Proxy n')) :: Int
    --           mklist i | i == 0 = pure []
    --                    | otherwise = (:) <$> xs <*> mklist (i-1)
    unfold xs = PLC1 . goodlength <$> unfold xs
        where n = fromInteger (natVal (Proxy :: Proxy n')) :: Int
              goodlength ys = assert (n > 0 && length ys == n) ys

instance (Applicative (ComonadicVector V.Vector), KnownNat n') =>
         Applicative (PLC1 n' c) where
  -- pure x = PLC1 (pure x)
  -- Note: 'pure' should use the neutral element of the monoid used
  -- for 'pos_' in '(<*>)'
  pure x = PLC1 (ComonadicVector (n - 1) (V.replicate n x))
      where n = fromInteger (natVal (Proxy :: Proxy n'))
  PLC1 fs <*> PLC1 xs = PLC1 (fs <*> xs)

-- instance (Alternative (ComonadicVector V.Vector), KnownNat n') =>
--          Alternative (PLC1 n' c) where
--   empty = PLC1 empty
--   PLC1 xs <|> PLC1 ys = PLC1 (xs <|> ys)

-- instance Monad (ComonadicVector V.Vector) => Monad (PLC1 n' c) where
--     -- (>>=) :: PLC1 n' c s -> (s -> PLC1 n' c t) -> PLC1 n' c t
--     PLC1 x >>= f = PLC1 (x >>= getPLC1 . f)
--
-- instance MonadPlus (ComonadicVector V.Vector) => MonadPlus (PLC1 n' c) where
--     mzero = PLC1 mzero
--     PLC1 xs `mplus` PLC1 ys = PLC1 (xs `mplus` ys)

-- instance Distributive (ComonadicVector V.Vector) =>
--          Distributive (PLC1 n' c) where
--   collect f arr = PLC1 (collect (getPLC1 . f) arr)

instance Comonad (ComonadicVector V.Vector) => Comonad (PLC1 n' c) where
  extract (PLC1 arr) = extract arr
  extend f (PLC1 arr) = PLC1 (extend (f . PLC1) arr)

instance ComonadStore Int (ComonadicVector V.Vector) =>
         ComonadStore Int (PLC1 n' c) where
  pos (PLC1 arr) = pos arr
  peek i (PLC1 arr) = peek i arr

instance (KnownNat n', QC.Arbitrary s) => QC.Arbitrary (PLC1 n' c s) where
  -- arbitrary = do
  --   arr <- QC.arbitrary
  --   return (PLC1 arr)
  -- shrink (PLC1 arr) = [PLC1 arr' | arr' <- QC.shrink arr]
  -- arbitrary = arbitraryDefault
  arbitrary = do
      xs <- V.generateM n (const QC.arbitrary)
      i <- QC.choose (0, V.length xs - 1)
      return (PLC1 (ComonadicVector i xs))
      where n = fromInteger (natVal (Proxy :: Proxy n'))
  shrink arr = [arr' | pos arr > 0] ++ traverse QC.shrink arr'
      where arr' = seek 0 arr

instance QC.CoArbitrary s => QC.CoArbitrary (PLC1 n' c s) where
    coarbitrary = QC.coarbitrary . getPLC1

instance QC.Function s => QC.Function (PLC1 n' c s) where
    function = QC.functionMap getPLC1 PLC1

-- instance Additive (ComonadicVector V.Vector) => Additive (PLC1 n' c) where
--     zero = PLC1 zero
--
-- instance (Additive (PLC1 n' c), Bounded c, RealFrac c) => Metric (PLC1 n' c) where
--     -- dot :: Num a => f a -> f a -> a
--     dot xs ys =
--         assert (length ys == n) $
--         case n of
--           0 -> 0
--           _ -> sum [w i * peek i xs * peek i ys | i <- [0 .. n-1]]
--         where n = length xs
--               xmin = minBound :: c
--               xmax = maxBound :: c
--               w i = realToFrac (weight n xmin xmax i)

instance ( KnownNat n'
         , AdditiveGroup (ComonadicVector V.Vector s)
         , AdditiveGroup s
         ) =>
         AdditiveGroup (PLC1 n' c s) where
  zeroV = pure zeroV
  negateV = fmap negateV
  (^+^) = liftA2 (^+^)

instance ( KnownNat n'
         , VectorSpace (ComonadicVector V.Vector s)
         , VectorSpace s
         ) =>
         VectorSpace (PLC1 n' c s) where
  type Scalar (PLC1 n' c s) = Scalar s
  (*^) x = fmap (x *^)

-- instance ( KnownNat n'
--          , Bounded c
--          , RealFrac c
--          , InnerSpace s
--          , VectorSpace (Scalar s)
--          , Fractional (Scalar (Scalar s))
--          ) =>
--          InnerSpace (PLC1 n' c s) where
--   -- xs <.> ys =
--   --   assert (length ys == length xs) $ sumV (weighted ((<.>) <$> xs <*> ys))

instance ( KnownNat n'
         , Integrable (PLC1 n')
         , IntegrableOk (PLC1 n') c s
         , Bounded c
         , RealFrac c
         , InnerSpace s
         , Fractional (Scalar s)
         ) => InnerSpace (PLC1 n' c s) where
  xs <.> ys = getPLC1 (extend weighted xs) <.> getPLC1 ys
    where
      weighted zs = weight n (xmin, xmax) (pos zs) *^ extract zs
      -- n = length xs
      n = fromInteger (natVal (Proxy :: Proxy n'))
      (xmin, xmax) = bounds xs

instance Function (PLC1 n') where
  type FunctionOk (PLC1 n') a b = ( KnownNat n'
                                  , Bounded a
                                  , RealFrac a
                                  , InnerSpace b
                                  , Floating (Scalar b)
                                  , Ord (Scalar b)
                                  , b ~ Double) -- TODO
  eval (PLC1 arr) x =
    case n of
      0 -> zeroV
      1 -> extract arr
      _ ->
        let y = ((x - cx) / rx + 1) / 2 * fromIntegral (n - 1)
            i = max 0 $ min (n - 2) $ floor y
            f0 = 1 - f1
            f1 = realToFrac (y - fromIntegral i)
        in f0 *^ peek i arr ^+^ f1 *^ peek (i + 1) arr
    where
      -- n = length arr
      n = fromInteger (natVal (Proxy :: Proxy n'))
      xmin = minBound `asTypeOf` x
      xmax = maxBound `asTypeOf` x
      cx = (xmax + xmin) / 2
      rx = (xmax - xmin) / 2

instance KnownNat n' => Discretization (PLC1 n') where
  discretized f = PLC1 (ComonadicVector 0 res)
    where
      n = fromInteger (natVal (Proxy :: Proxy n'))
      -- res = V.generate n go
      -- TODO: use (symmetric) triangular solver
      rhs = L.vector [go i | i <- [0..n-1]]
      mat = L.matrix n [ overlap n (xmin, xmax) i j
                       | i <- [0..n-1], j <- [0..n-1]]
      res = V.fromListN n $ L.toList $ mat L.<\> rhs
      (xmin, xmax) = (minBound, maxBound)
      _ = (f xmin, f xmax)      -- constrain types
      go i = integrate bf lo hi
        where
          bf x = b x *^ f x
          b x = basis n (xmin, xmax) i x
          (lo, hi) = support n (xmin, xmax) i

instance KnownNat n' => Integrable (PLC1 n') where
  type IntegrableOk (PLC1 n') a b = ()
  integral xs = sumV (extend weighted xs)
    where
      weighted ys = weight n (xmin, xmax) (pos ys) *^ extract ys
      -- n = length xs
      n = fromInteger (natVal (Proxy :: Proxy n'))
      (xmin, xmax) = bounds xs

instance Differentiable (PLC1 n') where
  type DifferentiableOk (PLC1 n') a b = ()
  type Dir (PLC1 n') = ()
  boundary () xs = extend bndry xs
    where
      bndry ys
        | n == 1     = zeroV *^ extract ys -- or zeroV
        | i == 0     = negateV (extract ys ^/ realToFrac dx)
        | i == n - 1 = extract ys ^/ realToFrac dx
        | otherwise  = zeroV *^ extract ys -- or zeroV
        where
          i = pos ys
      -- n = length xs
      n = fromInteger (natVal (Proxy :: Proxy n'))
      (xmin, xmax) = bounds xs
      rx = (xmax - xmin) / 2
      hdx = rx / fromIntegral (n - 1)
      dx = 2 * hdx
  derivative () xs = extend deriv xs
    where
      deriv ys
        | n == 1     = zeroV *^ extract ys -- or zeroV
        | i == 0     = (peek (i + 1) ys ^-^ extract ys) ^/ realToFrac dx
        | i == n - 1 = (extract ys ^-^ peek (i - 1) ys) ^/ realToFrac dx
        | otherwise  =
            (peek (i + 1) ys ^-^ peek (i - 1) ys) ^/ realToFrac (2 * dx)
        where
          i = pos ys
      -- n = length xs
      n = fromInteger (natVal (Proxy :: Proxy n'))
      (xmin, xmax) = bounds xs
      rx = (xmax - xmin) / 2
      hdx = rx / fromIntegral (n - 1)
      dx = 2 * hdx

triangle :: (Num a, Ord a) => a -> a
triangle x = max 0 (1 - abs x)

typed :: f a b -> a -> a
typed _ = id

bounds :: Bounded a => f a b -> (a, a)
bounds xs = (typed xs minBound, typed xs maxBound)

-- The basis functions are defined to have a maximum of '1'
basis :: (RealFrac a, Fractional b) => Int -> (a, a) -> Int -> a -> b
basis n (xmin, xmax) i x | n == 1    = 1
                         | otherwise = realToFrac (triangle y)
  where
    cx  = (xmin + xmax) / 2
    rx  = (xmax - xmin) / 2
    hdx = rx / fromIntegral (n - 1)
    xi  = cx + hdx * fromIntegral (2 * i + 1 - n)
    dx  = 2 * hdx
    y   = (x - xi) / dx

support :: (Fractional a, Ord a) => Int -> (a, a) -> Int -> (a, a)
support n (xmin, xmax) i | n == 1    = (xmin, xmax)
                         | otherwise = (max xmin (xi - dx), min xmax (xi + dx))
  where
    cx  = (xmin + xmax) / 2
    rx  = (xmax - xmin) / 2
    hdx = rx / fromIntegral (n - 1)
    xi  = cx + hdx * fromIntegral (2 * i + 1 - n)
    dx  = 2 * hdx

-- Volume
weight :: (RealFrac a, Fractional b) => Int -> (a, a) -> Int -> b
weight n (xmin, xmax) i | n == 1     = realToFrac (2 * rx)
                        | i == 0     = realToFrac hdx
                        | i == n - 1 = realToFrac hdx
                        | otherwise  = realToFrac (2 * hdx)
  where
    rx  = (xmax - xmin) / 2
    hdx = rx / fromIntegral (n - 1)

overlap :: (RealFrac a, Fractional b) => Int -> (a, a) -> Int -> Int -> b
overlap n (xmin, xmax) i j
    | n == 1                           = realToFrac (2 * rx)
    | i == j && (i == 0 || i == n - 1) = realToFrac (2 / 3 * hdx)
    | i == j                           = realToFrac (4 / 3 * hdx)
    | abs (i - j) == 1                 = realToFrac (1 / 3 * hdx)
    | otherwise                        = 0
  where
    rx  = (xmax - xmin) / 2
    hdx = rx / fromIntegral (n - 1)