text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeOperators       #-}
module Test.Grenade.Network where

import           Control.Monad                   (guard)
import           Control.Monad.ST                (runST)
import           Data.Constraint
import           Data.Proxy
import           Data.Singletons.Prelude.List
import           Data.Singletons.TypeLits
import           Data.Typeable
import qualified Data.Vector.Storable            as VS
import qualified Data.Vector.Storable.Mutable    as VS (write)
import           Hedgehog
import qualified Hedgehog.Gen                    as Gen
import           Hedgehog.Internal.Property      (failWith)
import           Hedgehog.Internal.Source

import           Grenade

#if MIN_VERSION_base(4,11,0)
import           GHC.TypeLits                    hiding (natVal)
#else
import           GHC.TypeLits
#endif
#if MIN_VERSION_base(4,9,0)
import           Data.Kind                       (Type)
#endif
import           Data.Singletons
import           Data.Singletons.Prelude.Num     ((%*))
import           GHC.TypeLits                    (someNatVal)

import           Test.Grenade.Layers.Convolution
import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix
import           Test.Hedgehog.TypeLits

import           Numeric.LinearAlgebra           (flatten)
import           Numeric.LinearAlgebra.Static    (extract, norm_Inf)
import           Unsafe.Coerce

data SomeNetwork :: Type where
    SomeNetwork :: ( SingI shapes, SingI (Head shapes), SingI (Last shapes), Show (Network layers shapes) ) => Network layers shapes -> SomeNetwork

instance Show SomeNetwork where
  show (SomeNetwork net) = show net

-- | Generate a random network of a random type
--
-- This is slightly insane for a few reasons. Everything must be wrapped up
-- in a SomeNetwork.
genNetwork :: Gen SomeNetwork
genNetwork =
  Gen.recursive Gen.choice [
      do SomeSing ( r :: Sing final ) <- genShape
         withSingI r $
           pure (SomeNetwork (NNil :: Network '[] '[ final ] ))
    ] [
      do SomeNetwork ( rest :: Network layers shapes ) <- genNetwork
         case ( sing :: Sing shapes ) of
           SNil -> Gen.discard -- Can't occur
           SCons ( h :: Sing h ) ( _ :: Sing hs ) ->
             withSingI h $
              case h of
                 D1Sing l@SNat ->
                   Gen.choice [
                     pure (SomeNetwork (Tanh    :~> rest :: Network ( Tanh    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Logit   :~> rest :: Network ( Logit   ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Relu    :~> rest :: Network ( Relu    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Elu     :~> rest :: Network ( Elu     ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Softmax :~> rest :: Network ( Softmax ': layers ) ( h ': h ': hs )))
                   , do -- Reshape to two dimensions
                        let divisors :: Integer -> [Integer]
                            divisors n = 1 : [x | x <- [2..(n-1)], n `rem` x == 0]
                        let len = fromIntegral $ natVal l
                        rs  <- Gen.element $ divisors len
                        let cs  = len `quot` rs
                        case ( someNatVal rs, someNatVal cs, someNatVal len ) of
                          ( Just (SomeNat (rs' :: Proxy inRows)), Just (SomeNat (cs' :: Proxy inCols)), Just (SomeNat (_ :: Proxy outLen ) )) ->
                            let p1 = singByProxy rs'
                                p2 = singByProxy cs'
                            in  case ( p1 %* p2, unsafeCoerce (Dict :: Dict ()) :: Dict ((inRows * inCols) ~ outLen), unsafeCoerce (Dict :: Dict ()) :: Dict (( 'D1 outLen ) ~ h )) of
                                  ( SNat, Dict, Dict ) ->
                                    pure (SomeNetwork (Reshape :~> rest :: Network ( Reshape ': layers ) ( ('D2 inRows inCols) ': h ': hs )))
                          _ -> Gen.discard -- Doesn't occur
                   ]
                 D2Sing r@SNat c@SNat ->
                   Gen.choice [
                     pure (SomeNetwork (Tanh    :~> rest :: Network ( Tanh    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Logit   :~> rest :: Network ( Logit   ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Relu    :~> rest :: Network ( Relu    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Elu     :~> rest :: Network ( Elu     ': layers ) ( h ': h ': hs )))
                   , do -- Build a convolution layer with one filter output
                        -- Figure out some kernel sizes which work for this layer
                        -- There must be a better way than this...
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c

                        let ok extent kernel = [stride | stride <- [ 1 .. extent ], (extent - kernel) `mod` stride == 0]

                        -- Get some kernels which will fit
                        kernel_r <- choose 1 output_r
                        kernel_c <- choose 1 output_c

                        -- Build up some strides which also fit
                        stride_r <- Gen.element $ ok output_r kernel_r
                        stride_c <- Gen.element $ ok output_c kernel_c

                        -- Determine the input size
                        let input_r = (output_r - 1) * stride_r + kernel_r
                        let input_c = (output_c - 1) * stride_c + kernel_c

                        guard (input_r < 100)
                        guard (input_c < 100)

                        -- Remake types for input
                        case ( someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal kernel_r, someNatVal kernel_c, someNatVal stride_r, someNatVal stride_c ) of
                          ( Just (SomeNat (_    :: Proxy inRows)),     Just (SomeNat  (_    :: Proxy inCols)),
                            Just (SomeNat (_    :: Proxy outRows)),    Just (SomeNat  (_    :: Proxy outCols)),
                            Just (SomeNat (pkr  :: Proxy kernelRows)), Just (SomeNat  (pkc  :: Proxy kernelCols)),
                            Just (SomeNat (_    :: Proxy strideRows)), Just (SomeNat  (_    :: Proxy strideCols))) ->
                              let p1 = singByProxy pkr
                                  p2 = singByProxy pkc
                              in  case ( p1 %* p2
                                        -- Fake it till you make it.
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D2 outRows outCols ) ~ h ))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outRows - 1) * strideRows) ~ (inRows - kernelRows)))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outCols - 1) * strideCols) ~ (inCols - kernelCols)))) of
                                    (SNat, Dict, Dict, Dict) -> do
                                        conv <- genConvolution
                                        pure (SomeNetwork (conv :~> rest :: Network ( Convolution 1 1 kernelRows kernelCols strideRows strideCols ': layers ) ( ('D2 inRows inCols) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a convolution layer with one filter output
                        -- Figure out some kernel sizes which work for this layer
                        -- There must be a better way than this...
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c

                        let ok extent kernel = [stride | stride <- [ 1 .. extent ], (extent - kernel) `mod` stride == 0]

                        -- Get some kernels which will fit
                        kernel_r <- choose 1 output_r
                        kernel_c <- choose 1 output_c
                        channels <- choose 1 10

                        -- Build up some strides which also fit
                        stride_r <- Gen.element $ ok output_r kernel_r
                        stride_c <- Gen.element $ ok output_c kernel_c

                        -- Determine the input size
                        let input_r = (output_r - 1) * stride_r + kernel_r
                        let input_c = (output_c - 1) * stride_c + kernel_c

                        guard (input_r < 100)
                        guard (input_c < 100)

                        -- Remake types for input
                        case ( someNatVal channels, someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal kernel_r, someNatVal kernel_c, someNatVal stride_r, someNatVal stride_c ) of
                          ( Just (SomeNat (chan :: Proxy channels)),
                            Just (SomeNat (pinr :: Proxy inRows)),     Just (SomeNat  (_    :: Proxy inCols)),
                            Just (SomeNat (_    :: Proxy outRows)),    Just (SomeNat  (_    :: Proxy outCols)),
                            Just (SomeNat (pkr  :: Proxy kernelRows)), Just (SomeNat  (pkc  :: Proxy kernelCols)),
                            Just (SomeNat (_    :: Proxy strideRows)), Just (SomeNat  (_    :: Proxy strideCols))) ->
                              let p1 = singByProxy pkr
                                  p2 = singByProxy pkc
                                  p3 = singByProxy chan
                              in  case ( p1 %* p2 %* p3
                                        , singByProxy pinr %* singByProxy chan
                                        -- Fake it till you make it.
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D2 outRows outCols ) ~ h ))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outRows - 1) * strideRows) ~ (inRows - kernelRows)))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outCols - 1) * strideCols) ~ (inCols - kernelCols)))) of
                                    (SNat, SNat, Dict, Dict, Dict) -> do
                                        conv <- genConvolution
                                        pure (SomeNetwork (conv :~> rest :: Network ( Convolution channels 1 kernelRows kernelCols strideRows strideCols ': layers ) ( ('D3 inRows inCols channels) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Pooling layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c

                        let ok extent kernel = [stride | stride <- [ 1 .. extent ], (extent - kernel) `mod` stride == 0]

                        -- Get some kernels which will fit
                        kernel_r <- choose 1 output_r
                        kernel_c <- choose 1 output_c

                        -- Build up some strides which also fit
                        stride_r <- Gen.element $ ok output_r kernel_r
                        stride_c <- Gen.element $ ok output_c kernel_c

                        -- Determine the input size
                        let input_r = (output_r - 1) * stride_r + kernel_r
                        let input_c = (output_c - 1) * stride_c + kernel_c

                        guard (input_r < 100)
                        guard (input_c < 100)

                        -- Remake types for input
                        case ( someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal kernel_r, someNatVal kernel_c, someNatVal stride_r, someNatVal stride_c ) of
                          ( Just (SomeNat (_    :: Proxy inRows)),     Just (SomeNat  (_    :: Proxy inCols)),
                            Just (SomeNat (_    :: Proxy outRows)),    Just (SomeNat  (_    :: Proxy outCols)),
                            Just (SomeNat (pkr  :: Proxy kernelRows)), Just (SomeNat  (pkc  :: Proxy kernelCols)),
                            Just (SomeNat (_    :: Proxy strideRows)), Just (SomeNat  (_    :: Proxy strideCols))) ->
                              let p1 = singByProxy pkr
                                  p2 = singByProxy pkc
                              in  case ( p1 %* p2
                                        -- Fake it till you make it.
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D2 outRows outCols ) ~ h ))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outRows - 1) * strideRows) ~ (inRows - kernelRows)))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outCols - 1) * strideCols) ~ (inCols - kernelCols)))) of
                                    (SNat, Dict, Dict, Dict) ->
                                        pure (SomeNetwork (Pooling :~> rest :: Network ( Pooling kernelRows kernelCols strideRows strideCols ': layers ) ( ('D2 inRows inCols) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Pad layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c

                        pad_left   <- choose 0 (output_r - 1)
                        pad_right  <- choose 0 (output_r - 1 - pad_left)
                        pad_top    <- choose 0 (output_c - 1)
                        pad_bottom <- choose 0 (output_c - 1 - pad_top)

                        let input_r = output_r - pad_top - pad_bottom
                        let input_c = output_c - pad_left - pad_right

                        guard (input_r > 1)
                        guard (input_c > 1)

                        -- Remake types for input
                        case ( someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal pad_left, someNatVal pad_right, someNatVal pad_top, someNatVal pad_bottom ) of
                          ( Just (SomeNat (_    :: Proxy inputRows)),  Just (SomeNat  (_    :: Proxy inputColumns)),
                            Just (SomeNat (_    :: Proxy outputRows)), Just (SomeNat  (_    :: Proxy outputColumns)),
                            Just (SomeNat (_    :: Proxy padLeft)),    Just (SomeNat  (_    :: Proxy padRight)),
                            Just (SomeNat (_    :: Proxy padTop)),     Just (SomeNat  (_    :: Proxy padBottom))) ->
                              case (  (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D2 outputRows outputColumns) ~ h ))
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inputRows + padTop + padBottom) ~ outputRows) )
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inputColumns + padLeft + padRight) ~ outputColumns ))) of
                                    ( Dict, Dict, Dict ) ->
                                        pure (SomeNetwork (Pad :~> rest :: Network ( Pad padLeft padTop padRight padBottom ': layers ) ( ('D2 inputRows inputColumns) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Crop layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c

                        crop_left   <- choose 0 10
                        crop_right  <- choose 0 10
                        crop_top    <- choose 0 10
                        crop_bottom <- choose 0 10

                        let input_r = output_r + crop_top + crop_bottom
                        let input_c = output_c + crop_left + crop_right

                        -- Remake types for input
                        case ( someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal crop_left, someNatVal crop_right, someNatVal crop_top, someNatVal crop_bottom ) of
                          ( Just (SomeNat (_    :: Proxy inputRows)),  Just (SomeNat  (_    :: Proxy inputColumns)),
                            Just (SomeNat (_    :: Proxy outputRows)), Just (SomeNat  (_    :: Proxy outputColumns)),
                            Just (SomeNat (_    :: Proxy cropLeft)),   Just (SomeNat  (_    :: Proxy cropRight)),
                            Just (SomeNat (_    :: Proxy cropTop)),    Just (SomeNat  (_    :: Proxy cropBottom))) ->
                              case (  (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D2 outputRows outputColumns) ~ h ))
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (outputRows + cropTop + cropBottom) ~ inputRows ) )
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (outputColumns + cropLeft + cropRight) ~ inputColumns ))) of
                                    (Dict, Dict, Dict) ->
                                        pure (SomeNetwork (Crop :~> rest :: Network ( Crop cropLeft cropTop cropRight cropBottom ': layers ) ( ('D2 inputRows inputColumns) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   ]
                 D3Sing r@SNat c@SNat f@SNat ->
                   Gen.choice [
                     pure (SomeNetwork (Tanh    :~> rest :: Network ( Tanh    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Logit   :~> rest :: Network ( Logit   ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Relu    :~> rest :: Network ( Relu    ': layers ) ( h ': h ': hs )))
                   , pure (SomeNetwork (Elu     :~> rest :: Network ( Elu     ': layers ) ( h ': h ': hs )))
                   , do -- Build a convolution layer with one filter output
                        -- Figure out some kernel sizes which work for this layer
                        -- There must be a better way than this...
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c
                        let output_f = fromIntegral $ natVal f

                        let ok extent kernel = [stride | stride <- [ 1 .. extent ], (extent - kernel) `mod` stride == 0]

                        -- Get some kernels which will fit
                        kernel_r <- choose 1 output_r
                        kernel_c <- choose 1 output_c
                        channels <- choose 1 10

                        -- Build up some strides which also fit
                        stride_r <- Gen.element $ ok output_r kernel_r
                        stride_c <- Gen.element $ ok output_c kernel_c

                        -- Determine the input size
                        let input_r = (output_r - 1) * stride_r + kernel_r
                        let input_c = (output_c - 1) * stride_c + kernel_c

                        guard (input_r < 100)
                        guard (input_c < 100)

                        -- Remake types for input
                        case ( someNatVal channels, someNatVal output_f, someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal kernel_r, someNatVal kernel_c, someNatVal stride_r, someNatVal stride_c ) of
                          ( Just (SomeNat (chan :: Proxy channels)),   Just (SomeNat  (_    :: Proxy filters)),
                            Just (SomeNat (pinr :: Proxy inRows)),     Just (SomeNat  (_    :: Proxy inCols)),
                            Just (SomeNat (_    :: Proxy outRows)),    Just (SomeNat  (_    :: Proxy outCols)),
                            Just (SomeNat (pkr  :: Proxy kernelRows)), Just (SomeNat  (pkc  :: Proxy kernelCols)),
                            Just (SomeNat (_    :: Proxy strideRows)), Just (SomeNat  (_    :: Proxy strideCols))) ->
                              let p1 = singByProxy pkr
                                  p2 = singByProxy pkc
                                  p3 = singByProxy chan
                              in  case ( p1 %* p2 %* p3
                                        , singByProxy pinr %* singByProxy chan
                                        -- Fake it till you make it.
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D3 outRows outCols filters) ~ h ))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outRows - 1) * strideRows) ~ (inRows - kernelRows)))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outCols - 1) * strideCols) ~ (inCols - kernelCols)))) of
                                    (SNat, SNat, Dict, Dict, Dict) -> do
                                        conv <- genConvolution
                                        pure (SomeNetwork (conv :~> rest :: Network ( Convolution channels filters kernelRows kernelCols strideRows strideCols ': layers ) ( ('D3 inRows inCols channels) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Pooling layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c
                        let output_f = fromIntegral $ natVal f

                        let ok extent kernel = [stride | stride <- [ 1 .. extent ], (extent - kernel) `mod` stride == 0]

                        -- Get some kernels which will fit
                        kernel_r <- choose 1 output_r
                        kernel_c <- choose 1 output_c

                        -- Build up some strides which also fit
                        stride_r <- Gen.element $ ok output_r kernel_r
                        stride_c <- Gen.element $ ok output_c kernel_c

                        -- Determine the input size
                        let input_r = (output_r - 1) * stride_r + kernel_r
                        let input_c = (output_c - 1) * stride_c + kernel_c

                        guard (input_r < 100)
                        guard (input_c < 100)

                        -- Remake types for input
                        case ( someNatVal output_f, someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal kernel_r, someNatVal kernel_c, someNatVal stride_r, someNatVal stride_c ) of
                          ( Just (SomeNat (chan :: Proxy filters)),
                            Just (SomeNat (pinr :: Proxy inRows)),     Just (SomeNat  (_    :: Proxy inCols)),
                            Just (SomeNat (_    :: Proxy outRows)),    Just (SomeNat  (_    :: Proxy outCols)),
                            Just (SomeNat (pkr  :: Proxy kernelRows)), Just (SomeNat  (pkc  :: Proxy kernelCols)),
                            Just (SomeNat (_    :: Proxy strideRows)), Just (SomeNat  (_    :: Proxy strideCols))) ->
                              let p1 = singByProxy pkr
                                  p2 = singByProxy pkc
                                  p3 = singByProxy chan
                              in  case ( p1 %* p2 %* p3
                                        , singByProxy pinr %* singByProxy chan
                                        -- Fake it till you make it.
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D3 outRows outCols filters) ~ h ))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outRows - 1) * strideRows) ~ (inRows - kernelRows)))
                                        , (unsafeCoerce (Dict :: Dict ()) :: Dict (((outCols - 1) * strideCols) ~ (inCols - kernelCols)))) of
                                    (SNat, SNat, Dict, Dict, Dict) ->
                                        pure (SomeNetwork (Pooling :~> rest :: Network ( Pooling kernelRows kernelCols strideRows strideCols ': layers ) ( ('D3 inRows inCols filters) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Pad layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c
                        let output_f = fromIntegral $ natVal f

                        pad_left   <- choose 0 (output_r - 1)
                        pad_right  <- choose 0 (output_r - 1 - pad_left)
                        pad_top    <- choose 0 (output_c - 1)
                        pad_bottom <- choose 0 (output_c - 1 - pad_top)

                        let input_r = output_r - pad_top - pad_bottom
                        let input_c = output_c - pad_left - pad_right

                        guard (input_r > 1)
                        guard (input_c > 1)

                        -- Remake types for input
                        case ( someNatVal output_f, someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal pad_left, someNatVal pad_right, someNatVal pad_top, someNatVal pad_bottom ) of
                          ( Just (SomeNat (chan :: Proxy filters)),
                            Just (SomeNat (pinr :: Proxy inputRows)),  Just (SomeNat  (_    :: Proxy inputColumns)),
                            Just (SomeNat (_    :: Proxy outputRows)), Just (SomeNat  (_    :: Proxy outputColumns)),
                            Just (SomeNat (_    :: Proxy padLeft)),    Just (SomeNat  (_    :: Proxy padRight)),
                            Just (SomeNat (_    :: Proxy padTop)),     Just (SomeNat  (_    :: Proxy padBottom))) ->
                              case (   singByProxy pinr %* singByProxy chan
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D3 outputRows outputColumns filters) ~ h ))
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inputRows + padTop + padBottom) ~ outputRows) )
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (inputColumns + padLeft + padRight) ~ outputColumns ))) of
                                    (SNat, Dict, Dict, Dict) ->
                                        pure (SomeNetwork (Pad :~> rest :: Network ( Pad padLeft padTop padRight padBottom ': layers ) ( ('D3 inputRows inputColumns filters) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                   , do -- Build a Crop layer
                        let output_r = fromIntegral $ natVal r
                        let output_c = fromIntegral $ natVal c
                        let output_f = fromIntegral $ natVal f

                        crop_left   <- choose 0 10
                        crop_right  <- choose 0 10
                        crop_top    <- choose 0 10
                        crop_bottom <- choose 0 10

                        let input_r = output_r + crop_top + crop_bottom
                        let input_c = output_c + crop_left + crop_right

                        -- Remake types for input
                        case ( someNatVal output_f, someNatVal input_r, someNatVal input_c, someNatVal output_r, someNatVal output_c, someNatVal crop_left, someNatVal crop_right, someNatVal crop_top, someNatVal crop_bottom ) of
                          ( Just (SomeNat (chan :: Proxy filters)),
                            Just (SomeNat (pinr :: Proxy inputRows)),  Just (SomeNat  (_    :: Proxy inputColumns)),
                            Just (SomeNat (_    :: Proxy outputRows)), Just (SomeNat  (_    :: Proxy outputColumns)),
                            Just (SomeNat (_    :: Proxy cropLeft)),   Just (SomeNat  (_    :: Proxy cropRight)),
                            Just (SomeNat (_    :: Proxy cropTop)),    Just (SomeNat  (_    :: Proxy cropBottom))) ->
                              case ( singByProxy pinr %* singByProxy chan
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( ( 'D3 outputRows outputColumns filters) ~ h ))
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (outputRows + cropTop + cropBottom) ~ inputRows ) )
                                    , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (outputColumns + cropLeft + cropRight) ~ inputColumns ))) of
                                    (SNat, Dict, Dict, Dict) ->
                                        pure (SomeNetwork (Crop :~> rest :: Network ( Crop cropLeft cropTop cropRight cropBottom ': layers ) ( ('D3 inputRows inputColumns filters) ': h ': hs )))
                          _ -> Gen.discard -- Can't occur
                     ]
    ]

-- | Test a partial derivative numerically for a random network and input
--
-- This is the most important test.
prop_auto_diff :: Property
prop_auto_diff = withDiscards 1000 . withTests 10000 . property $ do
  SomeNetwork (network :: Network layers shapes) <- forAll genNetwork
  (input  :: S (Head shapes))     <- forAllWith nice genOfShape
  (target :: S (Last shapes))     <- forAllWith nice oneUp
  (tested :: S (Head shapes))     <- forAllWith nice oneUp

  let (!tapes, !output) = runNetwork network input
  let (_, !backgrad)    = runGradient network tapes target
  let inputDiff         = input + tested * 0.00001
  let expected          = maxVal ( backgrad * tested )
  let (_, !outputDiff)  = runNetwork network inputDiff
  let result            = maxVal ( outputDiff * target - output * target ) / 0.00001

  result ~~~ expected

-- Make a shape where all are 0 except for 1 value, which is 1.
oneUp :: forall shape. ( SingI shape ) => Gen (S shape)
oneUp =
  case ( sing :: Sing shape ) of
    D1Sing SNat ->
      let x = 0 :: S ( shape )
      in  case x of
            ( S1D x' ) -> do
              let ex = extract x'
              let len = VS.length ex
              ix <- choose 0 (len - 1)
              let nx = runST $ do ex' <- VS.thaw ex
                                  VS.write ex' ix 1.0
                                  VS.freeze ex'
              maybe Gen.discard pure . fromStorable $ nx

    D2Sing SNat SNat ->
      let x = 0 :: S ( shape )
      in  case x of
            ( S2D x' ) -> do
              let ex = flatten ( extract x' )
              let len = VS.length ex
              ix <- choose 0 (len - 1)
              let nx = runST $ do ex' <- VS.thaw ex
                                  VS.write ex' ix 1
                                  VS.freeze ex'
              maybe Gen.discard pure . fromStorable $ nx

    D3Sing SNat SNat SNat ->
      let x = 0 :: S ( shape )
      in  case x of
            ( S3D x' ) -> do
              let ex = flatten ( extract x' )
              let len = VS.length ex
              ix <- choose 0 (len - 1)
              let nx = runST $ do ex' <- VS.thaw ex
                                  VS.write ex' ix 1
                                  VS.freeze ex'
              maybe Gen.discard pure . fromStorable $ nx

maxVal :: S shape -> RealNum
maxVal ( S1D x ) = norm_Inf x
maxVal ( S2D x ) = norm_Inf x
maxVal ( S3D x ) = norm_Inf x

(~~~) :: (Monad m, HasCallStack) => RealNum -> RealNum -> PropertyT m ()
(~~~) x y =
  if abs (x - y) < precision then
    success
  else
    withFrozenCallStack $
    failWith Nothing $ unlines [
        "━━━ Not Simliar ━━━"
        , show x
        , show y
        ]
  where precision | nameF == show (typeRep (Proxy :: Proxy Float)) = 2e-2
                  | otherwise = 2e-4
infix 4 ~~~

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE AllowAmbiguousTypes, ConstraintKinds #-}

module PCA.GaussianProcess
       ( GaussianProcess (..)
       , GPTrainingData (..)
       , PosteriorSample (..)
       , gpToPosteriorSample
       , kernelGP
       ) where

import Universum hiding (transpose, Vector)

import Control.Lens (makeLenses)

import Data.Array.Repa
import Data.Array.Repa.Repr.Unboxed (Unbox)
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)
import System.Random (Random, mkStdGen)

import PCA.Types (InputObservations(..), Matrix, Vector)
import PCA.Util

data GaussianProcess a = GaussianProcess
    { -- | we define the gaussian process by some kernel function
      _kernelGP :: Vector D a -> Vector D a -> Matrix D a
    }

data GPTrainingData a = GPTrainingData
    { _inputTrain  :: Vector D a  -- ^ an input training data
    , _outputTrain :: Vector D a  -- ^ an output training data
    }

makeLenses ''GaussianProcess
makeLenses ''GPTrainingData

newtype PosteriorSample a = PosteriorSample
    { unSample :: Matrix D a  -- ^ an posterior sample
    }

-- | The constraint kind required for getting a posterior sample by a given GP
type GPConstraint a =
  ( Field a
  , Random a
  , Unbox a
  , Floating a
  , Eq a
  )

-- | The Main GP function: get a posterior sample by some kernel function, input observations and the training data

gpToPosteriorSample
  :: GPConstraint a
  => InputObservations a        -- ^ input observations
  -> GaussianProcess a          -- ^ a kernel function
  -> GPTrainingData a           -- ^ a training data
  -> Int                        -- ^ the number of samples
  -> Maybe (PosteriorSample a)  -- ^ a posterior functional prior
gpToPosteriorSample (InputObservations observe@(ADelayed (Z :. len) _)) gP trainingData sampleNumber = do
  -- | The kernel applied to input test points (so-called K_ss)
  let covarianceMatrix = kernel observe observe

  -- | The kernel applied to input training points (so-called K)
  let trainingKernel = kernel inputTrain' inputTrain'

  -- | The Cholesky decomposition applied to kernel of training points (:), so-called L
  let cholK = cholSH $ trainingKernel +^
                (smap (* 0.00005) . identD . size . extent $ inputTrain')

  -- | a covariance between test points and input training points (so-called K_s)
  let testPointMean = kernel observe inputTrain'

  -- | (the roots of L * x = K_s)
  cholKSolve <- delay <$> linearSolveS cholK testPointMean

  -- | Here we solve  alinear system for output training points
  cholKSolveOut <- delay <$> linearSolveS cholK (transposeMatrix $ toMatrix outputTrain' len)

  -- | Here we compute a mean
  let mean = (transposeMatrix cholKSolveOut) `mulD` cholKSolve

  -- | A posterior
  let postF' = cholSH $
                     covarianceMatrix +^
                     ((smap (* 1.0e-6) (identD len)) -^
                     (transposeMatrix cholKSolve) `mulD` cholKSolve)

  -- | A posterior sample
  return $ (PosteriorSample $ mean +^ (functionalPrior postF' sampleNumber))
    where
       kernel = gP ^. kernelGP
       inputTrain' = trainingData ^. inputTrain
       outputTrain' = trainingData ^. outputTrain
       mulD m n = delay $ m `mulS` n
       functionalPrior matrix@(ADelayed (Z :. rows :. _) _) numberSample =
         delay $ matrix `mulS` randomCoeffs
         where
           randomCoeffs = randomMatrixD (mkStdGen (-4)) (rows, numberSample)

ext: hs

max_stars_repo_name: serokell/PCA

max_issues_repo_name: serokell/PCA

max_forks_repo_name: serokell/PCA


text:
{- |
Description :  DNA Markov models (substitution models or mutation models)
Copyright   :  (c) Dominik Schrempf 2017
License     :  GPLv3

Maintainer  :  dominik.schrempf@gmail.com
Stability   :  unstable
Portability :  non-portable (not tested)

This module provides rate matrices for the most common DNA substitution models
(they will be added according to need). They are also used as mutation models by
the boundary mutation model.

* Changelog

-}

module DNAModel
  ( DNAModelSpec(..)
  , DNAModel(..)
  , Nuc
  , getDNAModelInfoStr
  , rateMatrix
  , dnaModelSpecGetStateFreqVec
  , dnaModelScale
  ) where

import           Data.List             (intercalate)
import           Numeric.LinearAlgebra
import           RateMatrix

-- | The nucleotides.
data Nuc = A | C | G | T deriving (Eq, Show, Read, Ord, Bounded, Enum)

nNucs :: Int
nNucs = 1 + fromEnum (maxBound :: Nuc)

-- | The DNA model specifications. E.g., for the HKY model, we only have kappa.
data DNAModelSpec = JC |
                    -- | HKY model with transition to transversion ratio kappa.
                    HKY Double StationaryDist |
                    -- | GTR model with five params and state frequency vector.
                    GTR Double Double Double Double Double StationaryDist

instance Show DNAModelSpec where
  show JC        = "JC"
  show (HKY k f) = "HKY[" ++ show k ++ "]" ++ show f
  show (GTR a b c d e f) = "GTR[" ++ paramsStr ++ "]" ++ show f
    where paramsStr = intercalate "," (map show [a,b,c,d,e])

-- | Extract stationary frequency vector of a DNA model specification 'DNAModelSpec'.
dnaModelSpecGetStateFreqVec :: DNAModelSpec -> StationaryDist
dnaModelSpecGetStateFreqVec JC                = vector $ replicate nNucs 0.25
dnaModelSpecGetStateFreqVec (HKY _ f)         = f
dnaModelSpecGetStateFreqVec (GTR _ _ _ _ _ f) = f

-- | A rate matrix of a DNA models is called DNAModel.
data DNAModel = DNAModel { dnaRateMatrix :: RateMatrix
                         , dnaModelSpec  :: DNAModelSpec }

-- | The matrix of exchangeabilities; diagonal entries are not set.
exchangeabilityMatrix :: DNAModelSpec -> Matrix R
exchangeabilityMatrix JC = (4><4)
  [ 0.0, 1.0, 1.0, 1.0
  , 1.0, 0.0, 1.0, 1.0
  , 1.0, 1.0, 0.0, 1.0
  , 1.0, 1.0, 1.0, 0.0 ]
exchangeabilityMatrix (HKY k _) = (4><4)
  [ 0.0, 1.0,   k, 1.0
  , 1.0, 0.0, 1.0,   k
  ,   k, 1.0, 0.0, 1.0
  , 1.0,   k, 1.0, 0.0 ]
exchangeabilityMatrix (GTR a b c d e _) = (4><4)
  [ 0.0,   a,   b,   c
  ,   a, 0.0,   d,   e
  ,   b,   d, 0.0, 1.0
  ,   c,   e, 1.0, 0.0 ]
-- exchangeabilityMatrix _     = error "Model not yet supported."

-- | DNA model rate matrix normalized so that one mutation happens per unit time.
rateMatrix :: DNAModelSpec -> DNAModel
rateMatrix s = DNAModel rm s
  where
    f    = dnaModelSpecGetStateFreqVec s
    exch = exchangeabilityMatrix s
    rm   = normalizeRates f $ setDiagonal $ exch <> diag f
-- rateMatrix _ _       = error "Model not yet supported."

-- | Report a specific model.
getDNAModelInfoStr :: DNAModel -> String
getDNAModelInfoStr (DNAModel m s) =
  case s of
    JC        -> unlines
      [ "JC model with rate matrix"
      , show m ]
    (HKY k f) -> unlines
      [ "HKY model"
      , reportRateMatrix
      , reportExchangeabilityMatrix
      , reportStateFreqVec f
      , "And a kappa value of: " ++ show k ]
    (GTR a b c d e f) -> unlines
      [ "GTR model"
      , reportRateMatrix
      , reportExchangeabilityMatrix
      , reportStateFreqVec f
      , "And parameters: " ++ intercalate ", " (map show [a,b,c,d,e]) ]
      -- _ -> error "Model not yet supported."
  where reportRateMatrix   = "Rate matrix " ++ show m
        reportStateFreqVec f = "This corresponds to state frequencies (A, C, G, T): " ++ show f
        reportExchangeabilityMatrix = "Exchangeability matrix (diagonal set to 0.0) " ++
                                      show (exchangeabilityMatrix s)

-- | Scale a DNA model.
dnaModelScale :: Double -> DNAModel -> DNAModel
dnaModelScale s (DNAModel rm spec) = DNAModel (scale s rm) scaledSpec
  where scaledSpec = case spec of
                       JC              -> JC
                       HKY k f         -> HKY (s*k) f
                       GTR a b c d e f -> GTR (s*a) (s*b) (s*c) (s*d) (s*e) f

max_stars_repo_name: dschrempf/bmm-simulate

max_issues_repo_name: pomo-dev/bmm-simulate

max_forks_repo_name: pomo-dev/bmm-simulate


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TypeOperators #-}

-- |
-- Module      :  Spiral.Program.Monad
-- Copyright   :  (c) 2017 Drexel University
-- License     :  BSD-style
-- Maintainer  :  mainland@drexel.edu

module Spiral.Program.Monad (
    P(..),
    evalP,

    alwaysUnroll,
    shouldUnroll,

    comment,

    (.:=.),
    assignP,
    forP,
    tempP,

    cacheExp,

    newArray,
    cacheArray
  ) where

import Control.Monad.Exception (MonadException(..))
import Control.Monad.Fail (MonadFail)
import Control.Monad.IO.Class (MonadIO)
import Control.Monad.Primitive (PrimMonad(..))
import Control.Monad.Ref (MonadRef(..))
import Control.Monad.Reader (MonadReader(..),
                             ReaderT,
                             asks,
                             runReaderT)
import Control.Monad.State (MonadState(..),
                            StateT,
                            execStateT,
                            gets,
                            modify)
import Control.Monad.Trans.Class (MonadTrans(..))
import Data.Complex
import Data.IORef (IORef)
import Data.Map (Map)
import qualified Data.Map as Map
import Data.Ratio
import Data.Sequence ((|>))
import Text.PrettyPrint.Mainland

import Data.Heterogeneous
import Spiral.Array
import Spiral.Array.Repr.Hidden
import Spiral.Array.Repr.Complex
import Spiral.Array.Repr.Concrete
import Spiral.Config
import Spiral.Exp
import Spiral.Monad (MonadSpiral)
import Spiral.Program.Syntax
import Spiral.Util.Trace
import Spiral.Util.Uniq

data PEnv = PEnv { unroll :: Bool }

defaultPEnv :: PEnv
defaultPEnv = PEnv { unroll = False }

data PState = PState
    { decls  :: Decls
    , stms   :: Stms
    , ecache :: Map (Some Exp) Var
    }

defaultPState :: PState
defaultPState = PState
    { decls  = mempty
    , stms   = mempty
    , ecache = mempty
    }

newtype P m a = P { unP :: StateT PState (ReaderT PEnv m) a}
    deriving (Functor, Applicative, Monad, MonadFail, MonadIO,
              MonadException,
              MonadReader PEnv,
              MonadState PState,
              MonadUnique,
              MonadTrace,
              MonadConfig)

instance MonadTrans P where
    lift = P . lift . lift

deriving instance MonadRef IORef m => MonadRef IORef (P m)

instance PrimMonad m => PrimMonad (P m) where
    type PrimState (P m) = PrimState m
    primitive = P . primitive

instance MonadSpiral m => MonadSpiral (P m) where

evalP :: Monad m => P m () -> m Block
evalP m = do
    s <- runReaderT (execStateT (unP m) defaultPState) defaultPEnv
    return $ Block (decls s) (stms s)

appendDecl :: Monad m => Decl -> P m ()
appendDecl d = modify $ \s -> s { decls = decls s |> d }

appendStm :: Monad m => Stm -> P m ()
appendStm stm = modify $ \s -> s { stms = stms s |> stm }

inNewBlock :: Monad m => P m a -> P m (Block, a)
inNewBlock k = do
    old_decls <- gets decls
    old_stms  <- gets stms
    modify $ \s -> s { decls = mempty, stms = mempty }
    x <- k
    block <- Block <$> gets decls <*> gets stms
    modify $ \s -> s { decls = old_decls, stms = old_stms }
    return (block, x)

inNewBlock_ :: Monad m => P m a -> P m Block
inNewBlock_ m = fst <$> inNewBlock m

-- | Look up an expression in the expression cache.
lookupCachedExp :: Monad m => Exp a -> P m (Maybe Var)
lookupCachedExp e = Map.lookup (Some e) <$> gets ecache

-- | Insert an expression into the expression cache.
insertCachedExp :: Monad m => Var -> Exp a -> P m ()
insertCachedExp v e =
    modify $ \s -> s { ecache = Map.insert (Some e) v (ecache s) }

-- | Always unroll loops in the given computation.
alwaysUnroll :: Monad m => P m a -> P m a
alwaysUnroll = local $ \env -> env { unroll = True }

-- | Return 'True' if loops of the given size should be unrolled.
shouldUnroll :: MonadSpiral m => Int -> P m Bool
shouldUnroll n = do
    always <- asks unroll
    maxun  <- asksConfig maxUnroll
    return $ always || n <= maxun

-- | Comment on the code being generated.
comment :: MonadSpiral m => Doc -> P m ()
comment doc = whenDynFlag GenComments $ appendStm (CommentS doc)

-- | @'forP' begin end k@ generates code for a fro loop with bounds @begin@
-- and @end@, where the body of the loop is generated by calling theVarDVarD
-- continuation @k@ with an expression representing the loop variable.
forP :: MonadSpiral m => Int -> Int -> (Exp Int -> P m ()) -> P m ()
forP lo hi k = do
    should <- shouldUnroll (hi - lo)
    if should
      then mapM_ k [intE i | i <- [lo..hi-1::Int]]
      else do
        i     <- gensym "i"
        block <- inNewBlock_ $ k (VarE i)
        appendStm $ ForS i lo hi block

-- | Generate a temporary of the given type.
tempP :: forall a m . (Typed a, MonadSpiral m) => P m (Exp a)
tempP = do
    useComplexType <- asksConfig $ testDynFlag UseComplex
    go useComplexType tau
  where
    tau :: Type a
    tau = typeOf (undefined :: a)

    go :: Bool -> Type b -> P m (Exp b)
    go False (ComplexT tau') = do
        er <- go False tau'
        ei <- go False tau'
        return $ ComplexE er ei

    go _ tau = do
        v <- gensym "t"
        appendDecl $ VarD v tau
        return $ VarE v

-- | An alias for 'assignP'.
infix 4 .:=.
(.:=.) :: (Typed a, Num (Exp a), MonadSpiral m) => Exp a -> Exp a -> P m ()
(.:=.) = assignP

-- | Assign one expression to another.
assignP :: forall a m . (Typed a, Num (Exp a), MonadSpiral m) => Exp a -> Exp a -> P m ()
assignP e1 e2 = do
    e2' <- cacheExp e2
    assignS e1 e2'

assignS :: (Typed a, MonadSpiral m) => Exp a -> Exp a -> P m ()
assignS e1@(VarE v) e2 = do
    insertCachedExp v e2
    appendStm $ AssignS e1 e2

assignS (ComplexE e1r e1i) (ComplexE e2r e2i) = do
    assignS e1r e2r
    assignS e1i e2i

assignS e1 e2 =
    appendStm $ AssignS e1 e2

-- | Cache the given expression. This serves as a hint that it will be used
-- more than once.
cacheExp :: forall a m . (Typed a, Num (Exp a), MonadSpiral m) => Exp a -> P m (Exp a)
cacheExp e = do
    config <- askConfig
    cacheWithConfig config e

cacheWithConfig :: forall a m . (Typed a, Num (Exp a), MonadSpiral m)
                => Config
                -> Exp a
                -> P m (Exp a)
cacheWithConfig config = cache
  where
    useComplex, splitComplex, doCSE, doRewrite, difRewrite :: Bool
    useComplex   = UseComplex `testDynFlag` config
    splitComplex = SplitComplex `testDynFlag` config
    doCSE        = CSE `testDynFlag` config
    doRewrite    = Rewrite `testDynFlag` config
    difRewrite   = DifRewrite `testDynFlag` config

    cache :: forall a . (Typed a, Num (Exp a)) => Exp a -> P m (Exp a)
    cache e
      | doRewrite = rewrite e
      | otherwise = mustCache e

    rewrite :: forall a . (Typed a, Num (Exp a)) => Exp a -> P m (Exp a)
    rewrite (ConstE (FloatC x)) | x < 0 =
        return $ UnopE Neg (ConstE (FloatC (-x)))

    rewrite (ConstE (DoubleC x)) | x < 0 =
        return $ UnopE Neg (ConstE (DoubleC (-x)))

    rewrite e@ConstE{} =
        return e

    rewrite e@VarE{} =
        return e

    rewrite (ComplexE er ei) =
        ComplexE <$> rewrite er <*> rewrite ei

    -- Don't cache negated constants. We don't want to fall through to the next
    -- case because the call to 'negate' will re-negate the constant!
    rewrite e@(UnopE Neg ConstE{}) =
        return e

    -- Cache the term we are negating, not the negated term.
    rewrite (UnopE Neg e) =
        negate <$> rewrite e

    rewrite (UnopE op e) = do
        e' <- rewrite e
        mustCache (UnopE op e')

    rewrite (BinopE op e1 e2) = do
        e1' <- rewrite e1
        e2' <- rewrite e2
        go op e1' e2'
      where
        go :: Binop -> Exp a -> Exp a -> P m (Exp a)
        -- Choose canonical variable ordering
        go Add e1@VarE{} e2@VarE{} | e2 < e1 =
            mustCache (BinopE Add e2 e1)

        go Add e1@IdxE{} e2@IdxE{} | e2 < e1 =
            mustCache (BinopE Add e2 e1)

        go Sub e1@VarE{} e2@VarE{} | e2 < e1 =
            UnopE Neg <$> mustCache (BinopE Sub e2 e1)

        go Sub e1@IdxE{} e2@IdxE{} | e2 < e1 =
            UnopE Neg <$> mustCache (BinopE Sub e2 e1)

        go Mul e1@VarE{} e2@VarE{} | e2 < e1 =
            mustCache (BinopE Mul e2 e1)

        go Mul e1@IdxE{} e2@IdxE{} | e2 < e1 =
            mustCache (BinopE Mul e2 e1)

        -- Constants always come first
        go Mul e1 e2@ConstE{} | not (isConstE e1) =
            go Mul e2 e1

        -- Push negation out
        go Add e1 (UnopE Neg e2) =
            rewrite $ e1 - e2

        go Add (UnopE Neg e1) e2 =
            rewrite $ e2 - e1

        go Sub e1 (UnopE Neg e2) =
            rewrite $ e1 + e2

        go Sub (UnopE Neg e1) e2 =
            rewrite $ -(e1 + e2)

        go Mul e1 (UnopE Neg e2) =
            rewrite $ -(e1 * e2)

        go Mul (UnopE Neg e1) e2 =
            rewrite $ -(e1 * e2)

        -- Rewrite @w1*w2*t@ to @(w1*w2)*t@
        go Mul w1@(ConstE W{}) (BinopE Mul w2@(ConstE W{}) e) | difRewrite =
            cache $ (w1*w2) * e

        go Mul e1@(ConstE W{}) e2 =
            return $ BinopE Mul e1 e2

        -- Rewrite @w1*t1 + w2*t2@ to @w1*(t1 + (w2/w1)*t2)@ where either @w1@
        -- or @w2/w1@ is a 4th root of unity.
        go op e1@(BinopE Mul w1@(ConstE W{}) t1) e2@(BinopE Mul w2@(ConstE W{}) t2) | difRewrite, Just (+/-) <- fromAddSub op =
            case tau of
              ComplexT{} | isTrivialRoot w1 || isTrivialRoot (w2/w1) ->
                cache $ w1 * (t1 +/- ((w2/w1) * t2))
              _ -> mustCache (BinopE op e1 e2)
          where
            tau :: Type a
            tau = typeOf (undefined :: a)

        go op e1 e2 =
            mustCache (BinopE op e1 e2)

        fromAddSub :: Binop -> Maybe (Exp a -> Exp a -> Exp a)
        fromAddSub Add = return (+)
        fromAddSub Sub = return (-)
        fromAddSub _   = fail "isAddSub: neither (+) nor (-)"

    rewrite (ReE e) =
        ReE <$> rewrite e >>= mustCache

    rewrite (ImE e) =
        ImE <$> rewrite e >>= mustCache

    rewrite e =
        mustCache e

    -- When we are not using complex types, always separately cache the real and
    -- imaginary parts of a complex expression.
    mustCache :: forall a . (Typed a, Num (Exp a)) => Exp a -> P m (Exp a)
    mustCache (ComplexE er ei) =
        ComplexE <$> mustCache er <*> mustCache ei

    mustCache (ReE e) | not useComplex, splitComplex = do
        let (er, _ei) = forceUnComplexE e
        cache er

    mustCache (ImE e) | not useComplex, splitComplex = do
        let (_er, ei) = forceUnComplexE e
        cache ei

    mustCache e | not useComplex, splitComplex, ComplexT{} <- tau = do
        let (er, ei) = forceUnComplexE e
        ComplexE <$> cache er <*> cache ei >>= cache
      where
        tau :: Type a
        tau = typeOf (undefined :: a)

    mustCache e | doCSE = do
        maybe_v <- lookupCachedExp e
        case maybe_v of
          Just v  -> return $ VarE v
          Nothing -> do temp <- tempP
                        assignS temp e
                        return temp

    mustCache e =
        return e

-- | @isNthRoot n r@ returns 'True' if @r@ is an nth root of unity.
isNthRoot :: Int -> Exp (Complex a) -> Bool
isNthRoot n (ConstE (W n' k _)) = denominator ((k % n') * fromIntegral n) == 1
isNthRoot _ _                   = False

-- | 'isTrivialRoot' returns 'True' if its argument is a trivial root, i.e., a
-- fourth root of unity.
isTrivialRoot :: Exp (Complex a) -> Bool
isTrivialRoot = isNthRoot 4

-- | Create a new concrete array of the given shape.
newArray :: forall sh a m . (Shape sh, Typed a, Num (Exp a), MonadSpiral m)
         => (sh :. Int)
         -> P m (Array H (sh :. Int) (Exp a))
newArray (sh :. i) = do
    useComplexType <- asksConfig $ testDynFlag UseComplex
    v <- gensym "V"
    go useComplexType v (typeOf (undefined :: a))
  where
    go :: Bool -> Var -> Type a -> P m (Array H (sh :. Int) (Exp a))
    go False v (ComplexT tau) = do
        appendDecl $ ArrD v (sh :. 2*i) tau
        return $ H $ CMPLX $ C (sh :. 2*i) v

    go _ v tau = do
        appendDecl $ ArrD v (sh :. i) tau
        return $ H $ C (sh :. i) v

-- | Cache the contents of a matrix, returning it as a concrete matrix.
cacheArray :: forall a r sh m . (Typed a, Num (Exp a), Shape sh, IArray r (sh :. Int) (Exp a), MonadSpiral m)
           => Array r (sh :. Int) (Exp a)
           -> P m (Array H (sh :. Int) (Exp a))
cacheArray arr = do
    useComplexType <- asksConfig $ testDynFlag UseComplex
    v <- gensym "K"
    go useComplexType v (typeOf (undefined :: a))
  where
    go :: Bool -> Var -> Type a -> P m (Array H (sh :. Int) (Exp a))
    go False v ComplexT{} = do
        appendDecl $ ConstArrD v (RE arr)
        return $ H $ CMPLX $ C (extent (RE arr)) v

    go _ v _ = do
        appendDecl $ ConstArrD v arr
        return $ H $ C (extent arr) v

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
module Yasih.LispParser where

import Yasih.LispTypes

import Control.Monad
import Control.Monad.Except
import Data.Ratio
import Data.Complex
import Data.Array
import Data.IORef
import Numeric

import Text.Parsec.Char
import Text.ParserCombinators.Parsec

-- |Parse an expression
-- Parse and evaluate a LispVal returning a monadic value
readOrThrow :: Parser a -> String -> ThrowsError a
readOrThrow parser input =
    let checkEmpty = do
          rest <- getInput
          unless (null rest) $ fail $ "Non-empty: " ++ rest
        parser' = parser <* checkEmpty in
    case parse parser' "lisp" input of
        Left err -> throwError $ Parser err
        Right val -> return val

-- |Parse a single expression
readExpr :: String -> ThrowsError LispVal
readExpr = readOrThrow mainParser

-- |Parse multiple expressions
readExprList :: String -> ThrowsError [LispVal]
readExprList = readOrThrow $ many1 mainParser

-- | Discard leading whitespace
mainParser :: Parser LispVal
mainParser = do
    skipMany whitespace
    result <- parseExpr
    skipMany whitespace
    return result

-- |Parser that recognizes one of the symbols allowed in Scheme Ident.
symbol :: Parser Char
symbol = oneOf "!$%&|*/+-:<=?>@^_~"

whitespace :: Parser ()
whitespace =  void (void (try (space <|> tab)) <|> try parseComment)

tokenize :: Parser a -> Parser a
tokenize = (whitespace >>)

-- |Parser to ignore whitespace
spaces :: Parser ()
spaces = skipMany1 space

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

-- |Parse a boolean value
parseBool :: Parser LispVal
parseBool = do
    try $ char '#'
    (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

-- Exercise 2.2 and 2.3
-- Parse escaped Chars in strings
parseEscapedChars :: Parser Char
parseEscapedChars = do
    char '\\'
    x <- oneOf "\\\"nrt"
    return $ case x of
        '\\' -> x
        '"' -> x
        'n' -> '\n'
        'r' -> '\r'
        't' -> '\t'
        _ -> x

-- |Parse a string
parseString :: Parser LispVal
parseString = do
    char '"'
    x <- many (parseEscapedChars <|> noneOf "\"")
    char '"'
    return $ String x

-- |Parse a negative sign Char and return an
-- Integer value to be multiplied to a parsed number
parseIntegerSign :: Parser Integer
parseIntegerSign = do
    signChar <- optionMaybe $ oneOf "+-"
    return $ case signChar of
        Just '-' -> -1
        Just '+' -> 1
        _ -> 1

-- |Parse a negative sign Char and return an
-- Integer value to be multiplied to a parsed number
parseDoubleSign :: Parser Double
parseDoubleSign = do
    signChar <- optionMaybe $ oneOf "+-"
    return $ case signChar of
        Just '-' -> -1
        Just '+' -> 1
        _ -> 1

-- |Same as parseDoubleSign but fail if sign is not found
parseDoubleSignStrict :: Parser Double
parseDoubleSignStrict = do
    signChar <- optionMaybe $ oneOf "+-"
    case signChar of
        Just '-' -> return $ -1
        Just '+' -> return 1
        _ -> fail "no sign"

-- |Parse a Number
{-
    Parse many digits (Parser String)
    Then apply a LispVal Atom Number constructor
    Composed with read to the resulting value
    liftM is used to promote the function to a Monad
-}
parseNumber :: Parser LispVal
parseNumber = parseDecimal
    <|> parseDecimalExplBasis
    <|> parseBinary
    <|> parseOctal
    <|> parseHexadecimal

parseDecimal :: Parser LispVal
parseDecimal = do
    sign <- parseIntegerSign
    Number . (*sign) . read <$> many1 digit

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

bin2dig = bin2dig' 0
bin2dig' dig "" = dig
bin2dig' dig (x:xs) = bin2dig' (2 * dig + (if x == '0'
    then 0 else 1)) xs

parseBinary :: Parser LispVal
parseBinary = do
    try $ string "#b"
    sign <- parseIntegerSign
    x <- many1 $ oneOf "01"
    (return . Number . (*sign) . bin2dig) x

parseHexadecimal :: Parser LispVal
parseHexadecimal = do
    try $ string "#x"
    sign <- parseIntegerSign
    x <- many1 hexDigit
    (return . Number . (*sign) . fst . head . readHex) x

oct2dig  = ()
parseOctal :: Parser LispVal
parseOctal = do
    try $ string "#o"
    sign <- parseIntegerSign
    x <- many1 octDigit
    (return . Number . (*sign) . fst . head . readOct) x

parseFloat :: Parser LispVal
parseFloat = do
    sign <- parseDoubleSign
    x <- many1 digit
    char '.'
    y <- many1 digit
    (return . Float . (*sign) . fst . head . readFloat) $ x++"."++y

parseRatio :: Parser LispVal
parseRatio = do
    sign <- parseIntegerSign
    x <- many1 digit
    char '/'
    y <- many1 digit
    (return . Ratio) $ ((*sign) . read) x % read y

-- |Convert a LispVal Float or Integer to Haskell Double
toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromIntegral n

-- |Parse a Complex Number
parseComplex :: Parser LispVal
parseComplex = do
    realSign <- parseDoubleSign
    realVal <- try parseFloat <|> parseDecimal
    imagSign <- parseDoubleSignStrict
    imagVal <- try parseFloat <|> parseDecimal
    char 'i'
    (return . Complex) $ ((*realSign) . toDouble) realVal :+ ((*imagSign) . toDouble) imagVal

parseCharacter :: Parser LispVal
parseCharacter = do
    string "#\\"
    value <- try (string "newline" <|> string "space")
        <|> do
            x <- anyChar
            notFollowedBy alphaNum
            return [x]
    return $ Character $ case value of
        "space" -> ' '
        "newline" -> '\n'
        _ -> head value

-- |Parse an Expression
parseExpr :: Parser LispVal
parseExpr = try parseComplex
        <|> try parseRatio
        <|> try parseFloat
        <|> try parseNumber
        <|> try parseAtom
        <|> parseString
        <|> try parseBool
        <|> try parseCharacter
        <|> try parseQuoted
        <|> try parseQuasiQuoted
        <|> try parseUnQuote
        <|> try parseVector
        <|> try parseParens

-- |Parse a List of Atoms like a b c d
parseList :: Parser LispVal
parseList = List <$> sepEndBy parseExpr (skipMany1 whitespace)

-- |Parse a Dotted list (a b c . d)
parseDottedList :: Parser LispVal
parseDottedList = do
    -- Parse a List of 0 or more expressions
    -- Separated by whitespace
    head <- endBy parseExpr (skipMany1 whitespace)
    -- Then parse the remaining Expr after the dot
    tail <- char '.' >> skipMany whitespace >> parseExpr
    skipMany whitespace
    return $ DottedList head tail

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

-- |Parse a QuasiQuoted Expression
-- See https://schemers.org/Documents/Standards/R5RS/HTML/r5rs-Z-H-7.html#%_sec_4.2.6
parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted = do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

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

-- |Parse a list expression surrounded by parens
parseParens :: Parser LispVal
parseParens = do
    char '('
    skipMany space
    x <- try parseDottedList <|> try parseList
    char ')'
    return x

-- |Parse a Vector #(a b c)
-- |See https://schemers.org/Documents/Standards/R5RS/HTML/r5rs-Z-H-9.html#%_sec_6.3.6
parseVector :: Parser LispVal
parseVector = do
    string "#("
    skipMany space
    vals <- sepEndBy parseExpr whitespace
    char ')'
    return $
        Vector $ listArray (0, length vals -1) vals

-- |Parse a comment like
-- ; this is a lisp comment
parseComment :: Parser ()
parseComment =
    void (try (char ';') >> manyTill anyChar (void (char '\n') <|> eof))
    <|> void (try (string "#|") >> manyTill anyChar (void (string "|#")))


text:
module Day10.Day10 where

import           Data.Complex
import qualified Data.Set                      as S
import           Data.List
import           Lib
import qualified Data.List.Key                 as K
import           Data.Ord

type Asteroid = Complex Float

toCoords :: [String] -> [(Float, Float)]
toCoords input =
  [ (x, y) | (row, y) <- zip input [0 ..], (s, x) <- zip row [0 ..], s == '#' ]

parse :: [(Float, Float)] -> [Asteroid]
parse coords = [ x :+ y | (x, y) <- coords ]

genPolarCoords :: [Asteroid] -> [([Complex Float], Complex Float)]
genPolarCoords asteroids = map
  (\asteroid1 ->
    ( filter (\x -> x /= 0 :+ 0)
      $ map (\asteroid2 -> asteroid2 - asteroid1) asteroids
    , asteroid1
    )
  )
  asteroids

highest :: [([Complex Float], Complex Float)] -> ([Complex Float], Complex Float)
highest = K.maximum (length . S.fromList . map normalizePhase . fst)

normalizePhase :: Complex Float -> Float
normalizePhase x =
  let y = phase x
  in  if y < (-0.5) * pi
        then 2.5 * pi + y
        else 0.5 * pi + y

day10 :: IO ()
day10 = do
  contents <- fileContents "/src/Day10/input.txt"
  let coords               = toCoords (lines contents)
      asteroids            = genPolarCoords $ parse coords
      (bestStation, index) = highest asteroids
      sortOrder = mconcat [comparing normalizePhase, comparing magnitude]
      rows                 = length $ head $ lines contents
      groupedAsteroids =
        groupBy (\x y -> normalizePhase x == normalizePhase y)
          $ sortBy sortOrder bestStation
  print $ length groupedAsteroids
  print $ index + (head $ groupedAsteroids !! 199)


text:
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RankNTypes #-}
module Primitives.Math (primitives) where

import Data.Complex (Complex(..))
import Data.Functor ((<&>))
import Data.List (foldl', foldl1')
import Data.Ratio (approxRational)

import Val
import EvaluationMonad (throwError, panic)
import Types.Unwrappers (unwrapNum, unwrapRealNum, unwrapExactInteger)
import Primitives.Bool (predicate, predicateM)
import Primitives.Misc (valuesB)

addP, mulP, modP, quotP, remP :: Primitive
addP = numericBinop "+" (+) 0
mulP = numericBinop "*" (*) 1
modP = realBinop1 "mod" mod
quotP = realBinop1 "quotient" quot
remP = realBinop1 "remainder" rem

-- | Converts a (Haskell) builtin function on Integers to a 'Primitive'.
--   The primitive can only cause type errors;
--   Too many arguments will be folded over the operation;
--   none will use the default value.
numericBinop :: Symbol
             -> (Number -> Number -> Number)
             -> Number
             -> Primitive
numericBinop name op start = Prim name (AtLeast 0) $
  fmap (Number . foldl' op start) . mapM unwrapNum

-- | Same as numericBinop, except at least one arg must be given.
realBinop1 :: Symbol
              -> (RealNumber -> RealNumber -> RealNumber)
              -> Primitive
realBinop1 name op = Prim name (AtLeast 1) $
  fmap (Number . Real . foldl1' op) . mapM unwrapRealNum

guardDivZero :: Primitive -> Primitive
guardDivZero (Prim name arity f) = Prim name arity $ \args ->
    if any isExactZeroVal $ tail args
    then throwError $ Default "divide by zero"
    else f args
  where
    isExactZeroVal (Number n) = isExactZero n
    isExactZeroVal _ = False

-- | See the r7rs standard.
--
-- (- n) evaluates to -n, but (- x y z) evaluates to ((x - y) - z).
subP :: Primitive
subP = Prim "-" (AtLeast 1) $ \case
  [x]      -> Number . negate <$> unwrapNum x
  as@(_:_) -> Number . foldl1' (-) <$> mapM unwrapNum as
  _ -> panic "sub arity"

-- | r7rs: (/ n) is 1/n, (/ x y z) is ((x / y) / z)
divP :: Primitive
divP = guardDivZero $ Prim "/" (AtLeast 1) $ \case
  [] -> panic "div arity"
  [x] -> Number . (1 /) <$> unwrapNum x
  as  -> Number . foldl1' (/) <$> mapM unwrapNum as

negateP :: Primitive
negateP = Prim "negate" (Exactly 1) $
  \[x] -> Number . negate <$> unwrapNum x

numericPredicate :: Symbol -> (RealNumber -> Bool) -> Primitive
numericPredicate name p = predicateM name $ fmap p . unwrapRealNum

zeroCheck :: Primitive
zeroCheck = predicate "zero?" unwrapNum (== 0)

positiveCheck :: Primitive
positiveCheck = numericPredicate "positive?" (> 0)

negativeCheck :: Primitive
negativeCheck = numericPredicate "negative?" (< 0)

oddCheck :: Primitive
oddCheck = numericPredicate "odd?" $ (0 /=) . mod 2

evenCheck :: Primitive
evenCheck = numericPredicate "even?" $ (0 ==) . mod 2

-------------------------------------------------------------------------------
-- typically useful math operations
-------------------------------------------------------------------------------

inexactP :: Primitive
inexactP = Prim "inexact" (Exactly 1) $
  \[x] -> unwrapNum x <&> \case
    Real r -> Number $ Real $ realInexact r
    Complex (r :+ i) -> Number $ Complex $ realInexact r :+ realInexact i

rationalizeP :: Primitive
rationalizeP = Prim "rationalize" (Exactly 2) $
  \[x, eps] -> do
    rx <- unwrapRealNum x
    reps <- unwrapRealNum eps
    return $ Number $ Real $ Ratnum $ approxRational rx reps

squareP :: Primitive
squareP = Prim "square" (Exactly 1) $
  \[x] -> Number . (\n -> n*n) <$> unwrapNum x

-- | specialize ^ to Int powers
(^!) :: Num a => a -> Int -> a
x ^! n = x ^ n

exactIntegerSqrtP :: Primitive
exactIntegerSqrtP = Prim "exact-integer-sqrt" (Exactly 1) $
  \[k] -> do
    kN <- unwrapExactInteger k
    let root = squareRoot kN
        residual = kN - root ^! 2
    valuesB [makeBignum root, makeBignum residual]
  where
    -- This impl of the standard integer square root algorithm is from
    -- https://wiki.haskell.org/Generic_number_type#squareRoot
    squareRoot :: Integer -> Integer
    squareRoot 0 = 0
    squareRoot 1 = 1
    squareRoot n =
      let twopows = iterate (^! 2) 2
          (lowerRoot, lowerN) =
            last $ takeWhile ((n >=) . snd) $ zip (1:twopows) twopows
          newtonStep x = div (x + n `div` x) 2
          iters = iterate newtonStep (squareRoot (div n lowerN) * lowerRoot)
          isRoot r = r ^! 2 <= n && n < (r+1) ^! 2
       in head $ dropWhile (not . isRoot) iters

exptP :: Primitive
exptP = Prim "expt" (Exactly 2) $
  \[z1, z2] -> do
    z1' <- unwrapNum z1
    z2' <- unwrapNum z2
    return $ Number $ z1' ** z2'

-------------------------------------------------------------------------------
-- Complex functions
-------------------------------------------------------------------------------

makeRectangularP :: Primitive
makeRectangularP = Prim "make-rectangular" (Exactly 2) $
  \[x, y] -> do
    a <- unwrapRealNum x
    b <- unwrapRealNum y
    return $ Number $ Complex $ a :+ b

makePolarP :: Primitive
makePolarP = Prim "make-polar" (Exactly 2) $
  \[m, theta] -> do
    m' <- unwrapRealNum m
    theta' <- unwrapRealNum theta
    return $ Number $ Real m' * Complex (cos theta' :+ sin theta')

realPartP :: Primitive
realPartP = Prim "real-part" (Exactly 1) $
  \[z] -> unwrapNum z <&> \case
    Real r -> Number $ Real r
    Complex (r :+ _i) -> Number $ Real r

imagPartP :: Primitive
imagPartP = Prim "imag-part" (Exactly 1) $
  \[z] -> unwrapNum z <&> \case
    Real _r -> Number $ Real 0
    Complex (_r :+ i) -> Number $ Real i

magnitudeP :: Primitive
magnitudeP = Prim "magnitude" (Exactly 1) $
  \[z] -> Number . abs <$> unwrapNum z

angleP :: Primitive
angleP = Prim "angle" (Exactly 1) $
  \[z] -> unwrapNum z <&> \case
    Real r
      | r < 0 || isNegativeZero r -> Number pi
      | otherwise -> Number 0
    Complex (a :+ b) -> Number $ Real $ atan2 b a

liftFloatingMonop :: (Number -> Number) -> Builtin
liftFloatingMonop f = asBuiltin
  where asBuiltin [x] = do
          nx <- unwrapNum x
          return $ Number $ f nx
        asBuiltin _ = panic "liftFloatingMonop arity"

makeFloatingMonopPrim :: (Symbol, Number -> Number) -> Primitive
makeFloatingMonopPrim (name, f) = Prim name (Exactly 1) $ liftFloatingMonop f

floatingMonopPs :: [Primitive]
floatingMonopPs = map makeFloatingMonopPrim
  [ ("exp", exp)
  , ("sqrt", sqrt)
  , ("sin", sin)
  , ("cos", cos)
  , ("tan", tan)
  , ("asin", asin)
  , ("acos", acos)
  ]

logP :: Primitive
logP = Prim "log" (Between 1 2) $ \case
  [z] -> Number . log <$> unwrapNum z
  [z1, z2] -> do
    z1' <- unwrapNum z1
    z2' <- unwrapNum z2
    return $ Number $ logBase z2' z1'
  _ -> panic "logB arity"

atanP :: Primitive
atanP = Prim "atan" (Between 1 2) $ \case
  [z] -> Number . atan <$> unwrapNum z
  [y, x] -> do
    y' <- unwrapRealNum y
    x' <- unwrapRealNum x
    return $ Number $ Real $ atan2 y' x'
  _ -> panic "atanB arity"

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

primitives :: [Primitive]
primitives = [ addP
             , subP
             , mulP
             , divP
             , modP
             , quotP
             , remP
             , negateP
             , zeroCheck
             , positiveCheck
             , negativeCheck
             , oddCheck
             , evenCheck

             , inexactP
             , rationalizeP
             , atanP
             , exactIntegerSqrtP
             , squareP
             , exptP

             , makeRectangularP
             , makePolarP
             , realPartP, imagPartP
             , magnitudeP, angleP
             , logP
             ]
             ++ floatingMonopPs

max_stars_repo_name: JKTKops/Thcheme

max_issues_repo_name: JKTKops/Thcheme

max_forks_repo_name: JKTKops/Thcheme


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

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

import qualified Data.Text as T
import qualified Data.Text.IO as T
import qualified Data.Text.Read as T
import Data.Either (rights)
-------------------------
import Numeric.LinearAlgebra

readMMOM :: Int -> String -> IO [Double]
readMMOM nAtom foldernya = do
    fLLMF <- fmap (T.unpack . head) $ getLastLLMF foldernya
    putStrLn $ "===================processed LLMF=" ++ fLLMF
    mmom <- fmap (map T.double) $ inshell2text $ concat [ "mkdir -p temp; grep mmom ", fLLMF
                                            ,"| tail -n", show (nAtom + 1)
                                            ,"| head -n", show nAtom
                                            ,"| awk '{print $2}'"
                                          ]
    sdtMMOM <- fmap (map T.double) $ inshell2text $ concat [ "grep mmom ", fLLMF, "| grep ehf | tail -1 | sed -e 's/^.*mmom=//g'| awk '{print $1}'"
                                          ]
    return ( map fst $ rights $ concat [sdtMMOM,mmom])

----------------------------------------------------------------------
--getMMOM allArgs@(texFile:jd:jdHead:colAlign:xr:ymax':wTot:tumpuk:invS:tailer:foldernya:aos) = do

getMMOM :: [String] -> IO ()
getMMOM allArgs = do
--getMMOM allArgs = do
    putStrLn "===start ==== CPVO.IO.Reader.Ecalj.MMOM: getMMOM ==="
    --(invStat, ymax, xmin, xmax, ctrlAtoms, uniqAtoms, ctrlAtomicAOs,jdTable, cleanedJdHead, foldernya, tailer,colAlign,texFile) <- readHeaderData allArgs
    Right (invStat,_,_,_,ctrlAtoms,_, ctrlAtomicAOs,jdTable, cleanedJdHead,foldernya,tailer,_,texFile) <- readHeaderData allArgs

    -------------------------------generating data------------------------
    -------------------------------generating DOS data------------------------
    totalDOS <- readTotalDOSText tailer foldernya
    -------------------------------integrating DOS data-----------------------
    let intgTot = map (\i -> integrateToZero $ totalDOS ¿ [0,i]) $ flipBy invStat [1,2] -- run it on spin [1,2]
    putStrLn $ show intgTot
    -------------------------------integrating PDOS data------------------------
    let nAtom = length ctrlAtoms
    -------------------------------generating PDOS data------------------------
              -- map ditambah -1 karena input mengikuti gnuplot
              -- input : d kolom 6-10
              -- gnuplot : d kolom 6-10
              -- hmatrix : d kolom 5-9

  -------------------------------integrating PDOS data------------------------
    putStrLn $ "========invStat=" ++ (show invStat)
    (tMMomSD:mmomSD) <- fmap (map (* invStat)) $ readMMOM nAtom foldernya
    putStrLn $ show $ map (showDouble (3::Integer)) mmomSD
    putStrLn $ show tMMomSD
    putStrLn "==========show tMMomSD==========="
    pdosAtomicAll <- readPDOS invStat tailer foldernya ctrlAtomicAOs
    let integratedAtomicPDOS = integrateAtomicPDOS pdosAtomicAll
    let rIntgAll' = rendertable
         $ (:) cleanedJdHead
         $ (:) (concat [ ["Total" ]
           , ["  "]
           , map (showDouble (3::Integer)) $ (\[t,iu,idn] -> [t,iu-idn,t-(iu-idn)]) $ (tMMomSD:intgTot)
           ])
         $ zipWith (\a b -> a:b) (map show ([1,2..]::[Integer]))
         $ zipWith (\sdMom (intMom,(_,(j,_))) -> j:(map (showDouble (3::Integer)) [sdMom,intMom,sdMom-intMom])) mmomSD
         $ map (\(iu,idn,b) -> ((iu-idn),b) ) integratedAtomicPDOS
    let rIntgAll = unlines  [
                            rIntgAll'
                            , jdTable
                            ]
    {-
    resIntAll' <- markdownToTex rIntgAll
    let resIntAll = T.replace "\\}" "}"
                  $ T.replace "\\{" "{" $ T.pack
                  $ unlines [
                            "\\begin{longtable}[]{" ++ colAlign ++ "}"
                            , unlines $ tail $ lines $ T.unpack resIntAll'
                            ]
    putStrLn rIntgAll
    T.putStrLn resIntAll
    -}
    T.writeFile texFile $ T.pack rIntgAll
    putStrLn "===done CPVO.IO.Reader.Ecalj.MMOM: getMMOM ==="

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


text:
module Statistics.Distribution.LogNormal ( logNormalDistr
                                         , logNormalDistr'
                                         , getNormal
                                         ) where

import Statistics.Distribution
import Statistics.Distribution.Normal

data LogNormalDistribution = LND { getNormal :: NormalDistribution }
    deriving (Eq, Read, Show)

logNormalDistr :: Double -> Double -> LogNormalDistribution
logNormalDistr m s = LND $ normalDistr m s

-- gives the approximate mean and stdDev for the LogNormal.
logNormalDistr' :: Double -> Double -> LogNormalDistribution
logNormalDistr' m s = LND $ normalDistr mu sigma
    where mu = log $ m*m / sqrt (s*s + m*m)
          sigma = sqrt . log $ 1 + ((s*s) / (m*m))

instance Distribution LogNormalDistribution where
    cumulative (LND nd) = cumulative nd . log

instance Mean LogNormalDistribution where
    mean (LND nd) = exp $ variance nd / 2 + mean nd

instance MaybeMean LogNormalDistribution where
    maybeMean = Just . mean

instance Variance LogNormalDistribution where
    variance lnd@(LND nd) = (exp (variance nd) - 1) * (mean lnd ** 2)

instance MaybeVariance LogNormalDistribution where
    maybeVariance = Just . variance

instance ContGen LogNormalDistribution where
    genContVar (LND nd) g = exp <$> genContVar nd g

instance ContDistr LogNormalDistribution where
    quantile (LND nd) = quantile nd . log
    density (LND nd) x | x <= 0     = 0.0
                       | otherwise = density nd (log x) / (x * sqrt (stdDev nd) )

max_stars_repo_name: cspollard/HMCMC

max_issues_repo_name: cspollard/HMCMC

max_forks_repo_name: cspollard/HMCMC


text:
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedLists       #-}
{-# LANGUAGE ViewPatterns          #-}

-- | Main module
module SatOpt (main) where
--module SatOpt (module SatOpt) where

import           SatOpt.FFT                    as SatOpt
import           SatOpt.Optimize               as SatOpt
import           SatOpt.Render                 as SatOpt
import           SatOpt.Utility                as SatOpt
-- GENERATE: import New.Module as SatOpt

import           Data.Colour                   (blend)
import           Data.Colour.RGBSpace
import           Data.Colour.RGBSpace.HSL
import           Data.Complex
import           Data.List
import           Data.Packed.Matrix            (mapMatrix)
import           Debug.Trace                   (trace)
import           Numeric.LinearAlgebra.HMatrix hiding (disp, matrix, size)
import qualified Numeric.LinearAlgebra.HMatrix as HM

data Field2D a = FMat  (Matrix a)
               | FFT   (Matrix a)
               | FDisp (Matrix a)
               | FFunc (ℝ -> ℝ -> a)

data HSource = HPointSources  [(ℝ, ℝ, ℂ)]
             | HCircleSources ℝ [ℂ]
             | HSum [HSource]

hsrcToField :: HSource -> Field2D ℂ
hsrcToField (HPointSources ps)   = FFunc $ pointSources ps
hsrcToField (HCircleSources r p) = FFunc $ circleOfSources r p
hsrcToField (HSum [h])           = hsrcToField h
hsrcToField (HSum (h:hs))        = fieldSum (hsrcToField h) (hsrcToField (HSum hs))

fieldSum :: Field2D ℂ -> Field2D ℂ -> Field2D ℂ
fieldSum (FMat m1)  (FMat m2)  = FMat (m1 + m2)
fieldSum (FFunc f1) (FFunc f2) = FFunc (\x y -> f1 x y + f2 x y)
fieldSum (FMat m1)  (FFunc f2) = fieldSum (FFunc f2) (FMat m1)
fieldSum (FFunc f1) (FMat m2)  = fieldSum (matField (FFunc f1)) (FMat m2)

data HDisp = HConstDisp ℂ
           | HFuncDisp (ℝ -> ℝ -> ℂ)

hdspToField :: HDisp -> Field2D ℂ
hdspToField (HConstDisp p) = constF p
hdspToField (HFuncDisp f) = FFunc f

data HelmConf = HelmConf { helmsrc :: HSource
                         , helmdsp :: HDisp
                         }

data State a = State { amplitude  :: Field2D ℂ
                     , dispersion :: Field2D ℂ
                     , sources    :: Field2D ℂ
                     , time       :: ℝ
                     , freq       :: ℝ
                     , conf       :: a
                     , configure  :: State a -> State a
                     , concrete   :: Bool
                     }

constF :: Element a => a -> Field2D a
constF x = matField (FFunc (\_ _ -> x))

defState :: State a
defState = State { amplitude  = constF (0 :+ 0)
                 , dispersion = constF (1 :+ 0.1)
                 , sources    = constF (0 :+ 0)
                 , time       = 0
                 , freq       = 1
                 , conf       = undefined
                 , configure  = const defState
                 , concrete   = False
                 }

fromFieldC :: Field2D ℂ -> Matrix ℂ
fromFieldC (FFunc f) = sample2DM f [0, 1 .. size - 1] [0, 1 .. size - 1]
fromFieldC (FMat m)  = m
fromFieldC (FFT m)   = ifft2dM m
fromFieldC (FDisp m) = sqrt $ (recip m) - (square normm)

matField :: Element a => Field2D a -> Field2D a
matField (FFunc f) = FMat $ sample2DM f [0, 1 .. size - 1] [0, 1 .. size - 1]
matField x         = x

fftMat :: Field2D ℂ -> Matrix ℂ
fftMat (FFT m)  = m
fftMat (FMat m) = fft2dM m
fftMat f        = fftMat $ matField f

dispMat :: Field2D ℂ -> Matrix ℂ
dispMat (FDisp m) = m
dispMat (FMat m)  = recip $ square m + square normm
dispMat f         = dispMat $ matField f

makeConcrete :: (State a) -> (State a)
makeConcrete st@(State { concrete = True }) = st
makeConcrete (st@(State { dispersion = d, sources = s, time = t, freq = f }))
  = st { amplitude = ampl
       , dispersion = FDisp disp
       , sources = FFT srcs
       , concrete = True }
  where
    disp = dispMat $ matField d
    srcs = fftMat $ matField s
    ampl = FMat $ timeShift f t $ ifft2dM $ srcs * disp

phaseShift :: Matrix ℝ -> Matrix ℂ -> Matrix ℂ
phaseShift p m = mags * zipMatrixWith (:+) (cos angles) (sin angles)
  where
    angles = p + (mapMatrix realPart m `arctan2` mapMatrix imagPart m)
    mags = m * conj m

timeShift :: ℝ -> ℝ -> Matrix ℂ -> Matrix ℂ
timeShift f t m = case matField $ constF (tau * f * t) of
                   FMat p -> phaseShift p m
                   _      -> error "matField output invalid"

size :: Num a => a
size = 256

tau :: ℝ
tau = 2 * pi

rootsOfUnity :: Int -> [ℂ]
rootsOfUnity i = take i $ iterate (rotate (tau / fromIntegral i)) (1 :+ 0)

pointsOnCircle :: Int -> ℝ -> [(ℝ, ℝ)]
pointsOnCircle i r = map (fromC . scaleC r) $ rootsOfUnity i
  where
    fromC (a :+ b) = (a, b)
    scaleC s c = (s :+ 0) * c

-- | Take a position with respect to the center and move its origin to the top right
center :: Floating a => a -> a
center x = x - (size / 2)

pointSources :: [(ℝ, ℝ, ℂ)] -> ℝ -> ℝ -> ℂ
pointSources ps x y = sum $ map (\(cx, cy, s) -> s * delta cx cy x y) ps

circleOfSources :: ℝ -> [ℂ] -> ℝ -> ℝ -> ℂ
circleOfSources rad cs = pointSources $ zipCs $ pointsOnCircle (length cs) rad
  where
    zipCs = map (\(z, (x, y)) -> (x, y, z)) . zip cs

--const $ FFunc testFunc

-- testFunc :: HelmConf -> Field2D ℂ
-- testFunc (HelmConf (HPointSources ps) _)   = FFunc $ pointSources ps
-- testFunc (HelmConf (HCircleSources r p) _) = FFunc $ circleOfSources r p
  -- where
  --   pol m p = mkPolar m ((pi/180)*p)
-- testFunc = pointSources [ (0.48, 0.48, 0.5 :+ 0)
--                         , (0.52, 0.52, (-0.5) :+ 0)
--                         ]

-- 0.1 [ pol 1 0
--                                , pol 1 60
--                                , pol 1 120
--                                , pol 1 120
--                                , pol 1 60
--                                , pol 1 0
--                                ]

normm :: Matrix ℂ
--im, jm :: Matrix ℂ
--im = sample2DM (\x _ -> x :+ 0) [0, 1 .. size - 1] [0, 1 .. size - 1]
--jm = sample2DM (\_ y -> y :+ 0) [0, 1 .. size - 1] [0, 1 .. size - 1]
normm = sample2DM (\x y -> norm x y :+ 0) [0, 1 .. size - 1] [0, 1 .. size - 1]

eps :: ℝ
eps = 5

delta :: ℝ -> ℝ -> ℝ -> ℝ -> ℂ
delta ((*size) -> cx) ((*size) -> cy) ((*size) -> x) ((*size) -> y)
  | norm (center x - cx) (center y - cy) < size * eps    = 1
  | otherwise                                                       = 0

rendCmp_ :: ℂ -> (ℝ, ℝ, ℝ)
rendCmp_ (polar -> (m, p)) = uncurryRGB (\x y z -> (x, y, z)) $ hsl a 1 m
  where
--    a = if p < 0 then 90 else 0
    a = (p + pi) * (180 / pi)

rendCmp :: Matrix ℂ -> (Matrix ℝ, Matrix ℝ, Matrix ℝ)
rendCmp = matrify . unzip3 . map (unzip3 . map rendCmp_) . toLists
  where
    matrify (x, y, z) = (fromLists x, fromLists y, fromLists z)

rotate :: ℝ -> ℂ -> ℂ
rotate r = uncurry mkPolar . (\(m, p) -> (m, p + r)) . polar

visual :: (ℝ, ℝ, ℝ, ℝ) -> ℂ -> ℂ
visual (xa, xb, ya, yb) (polar -> (m, p))
  = mkPolar (((yb - ya) * (m - xa)/(xb - xa)) + ya) p

visualM :: Matrix ℂ -> Matrix ℂ
visualM cs@(fromComplex -> (r, i)) = mapMatrix (visual (small, big, 0, 1)) cs
  where
    !mags  = flatten $ sqrt $ (r * r) + (i * i)
    !small = minElement mags
    !big   = maxElement mags

updateState :: State a -> State a
--updateState = id
updateState s@State { time = t } = s { time = t + 0.01, concrete = False }

data RenderMode = RendAmp | RendSrc | RendDsp | RendAmpMag | RendAmpPhs

renderState :: RenderMode -> State a -> RGBTrips
renderState r = case r of
  RendAmp    -> rend . concAmp
  RendAmpMag -> rend . concAmpMag
  RendAmpPhs -> rend . concAmpPhs
  RendSrc    -> rend . concSrc
  RendDsp    -> rend . concDsp
  where
  wrap f = fromFieldC . f . makeConcrete
  rend =  rendCmp . visualM
  concAmp    = wrap amplitude
  concAmpMag = mag . fromComplex . concAmp
  concAmpPhs = phs . fromComplex . concAmp
  concSrc    = wrap sources
  concDsp    = wrap dispersion
  mag (x, y) = toComplex $ ((x*x) + (y*y), y - y)
  phs (x, y) = toComplex $ ((y `arctan2` x), y - y)

changeRenderMode :: RenderMode -> Conf (State a) -> Conf (State a)
changeRenderMode r c = c { render = renderState r }

getConf :: Conf (State HelmConf) -> HelmConf
getConf (Conf { state = (State { conf = cn })}) = cn

changeConf :: (HelmConf -> HelmConf) -> Conf (State HelmConf) -> Conf (State HelmConf)
changeConf h c@(Conf { state = s@(State { configure = f, conf = cn })})
  = c { state = configure s' s' }
  where
    s' = s { conf = h cn }

relPhase (polar -> (_, p1)) (polar -> (_, p2)) = abs $ p1 - p2

--getPhase c@(HelmConf { helmsrc = HCircleSources _ [c1, c2] }) = Just $ relPhase c1 c2
--getPhase _ = Nothing
getPhase c@(HelmConf { helmsrc = HPointSources [_, (x, _, _)] }) = Just x
getPhase _ = Nothing

changePhase :: Double -> Conf (State HelmConf) -> Conf (State HelmConf)
changePhase p = changeConf shift
  where
    shift c@(getPhase -> m) = case m of
                               Just ph -> c { helmsrc = testSource $ ph + p }
                               Nothing -> c

keyState :: KMState -> Conf (State HelmConf) -> IO (Conf (State HelmConf))
keyState ""  c = return c
keyState "d" c = return $ changeRenderMode RendDsp c
keyState "s" c = return $ changeRenderMode RendSrc c
keyState "a" c = return $ changeRenderMode RendAmp c
keyState "m" c = return $ changeRenderMode RendAmpMag c
keyState "p" c = return $ changeRenderMode RendAmpPhs c
keyState "w" c = do
  let p = changePhase 0.5 c
  putStrLn $ "Shifted by 0.5 to " ++ show (getPhase $ getConf c)
  return p
keyState "q" c = do
  let p = changePhase (-0.5) c
  putStrLn $ "Shifted by -0.5 to " ++ show (getPhase $ getConf c)
  return p
keyState x   c = return $ trace x c

pol m p = mkPolar m (p * pi / 180)

testSource :: Double -> HSource
testSource p = HPointSources [ (size/2, size/2, (-1) :+ 0)
                             , (p, p, 1 :+ 0)
                             ]
-- testSource p = HSum [ HCircleSources 0.2  [ pol 1 0
--                                           , pol 1 0
--                                           , pol 1 0
--                                           , pol 1 0
--                                           ]
--                     , HCircleSources 0.15 [ pol 1 0
--                                           , pol 1 90
--                                           , pol 1 180
--                                           , pol 1 270
--                                           ]
--                     ]

dspfun :: ℝ -> ℝ -> ℂ
dspfun x y
  | outbox (center x) (center y) = 1.2 :+ 0
  | otherwise                    = 1 :+ 0
  where
    test = 1
    outbox a b
      | a < (-32) = True
      | a > 32    = True
      | b < (-32) = True
      | b > 32    = True
      | otherwise = False

testConf :: Double -> HelmConf
testConf p = HelmConf { helmsrc = testSource p
                      , helmdsp = HFuncDisp dspfun
                      }

-- testConf = HelmConf { helmsrc = HCircleSources 0.1 [ pol 1 0
--                                                    , pol 1 0
--                                                    , pol 1 0
--                                                    , pol 1 0
--                                                    , pol 1 0
--                                                    , pol 1 0 ]
--                     , helmdsp = HConstDisp 0.001
--                     }
--   where
--     pol m p = mkPolar m (p * pi / 180)

testConfigure st@(State { conf = HelmConf { helmdsp = d, helmsrc = s }})
  = st { sources = hsrcToField s
       , dispersion = hdspToField d
       , concrete = False
       }

initState :: State HelmConf
initState = defState { freq = 3
                     , conf = testConf 0
                     , configure = testConfigure
                     }

config :: Conf (State HelmConf)
config = changeConf id
         Conf { keyBinds = keyState
              , state    = initState
              , evolve   = updateState
              , render   = renderState RendAmp
              , canvas   = (size, size)
              }

-- | TODO
main :: IO ()
--main = optimize
main = runDisp config

max_stars_repo_name: taktoa/SatOpt

max_issues_repo_name: taktoa/SatOpt

max_forks_repo_name: taktoa/SatOpt


text:
module Foreign.Storable.OrphansSpec (main, spec) where

import Test.Hspec
import Data.Complex
import Data.Orphans ()
import Data.Ratio
import Foreign.Storable

main :: IO ()
main = hspec spec

spec :: Spec
spec = do
  describe "Storable Complex instance" $ do
    it "has twice the sizeOf its realPart" $ do
      sizeOf (undefined :: Complex Double) `shouldBe` 2*sizeOf (1 :: Double)
    it "has the alignment of its realPart" $ do
      alignment (undefined :: Complex Double) `shouldBe` alignment (1 :: Double)

  describe "Storable Ratio instance" $ do
    it "has twice the sizeOf its parameterized type" $ do
      sizeOf (undefined :: Ratio Int) `shouldBe` 2*sizeOf (1 :: Int)
    it "has the alignment of its parameterized type" $ do
      alignment (undefined :: Ratio Int) `shouldBe` alignment (1 :: Int)

max_stars_repo_name: haskell-compat/base-orphans

max_issues_repo_name: haskell-compat/base-orphans

max_forks_repo_name: haskell-compat/base-orphans


text:
module Utils where

import Data.Vector.Unboxed as U
import Prelude as P
import Data.Decimal
import GHC.Word

import qualified Statistics.Matrix as M
import qualified Statistics.Matrix.Types as T
import qualified Statistics.Matrix.Algorithms as Alg
import qualified Statistics.Matrix.Function as Func
import qualified Statistics.Matrix.Mutable as Mut

import qualified Fixtures as F

-- | Operations on vectors

-- | Operations on matrices

-- | Check if two matrices are equal within some precision
isEqual :: Word8 -> T.Matrix -> T.Matrix -> Bool
isEqual prec a b =
    let
        a' = makeDecimal prec a
        b' = makeDecimal prec b
    in
        a' == b'

------------------------------ matrix property checks --------------------------------------

-- | Check if given matrix is orthogonal
isOrtho :: T.Matrix -> Bool
isOrtho mat = isEqual 8 (M.multiply mat (M.transpose mat)) F.matId

-- | Checks if the n'th row of the matrix have 'n' leading zeros
hasLeadingZeros :: Int -> T.Vector -> Bool
hasLeadingZeros n row = U.foldl (&&) True (U.map (== 0) (U.take n row))

-- | pass the n'th row of the matrix to hasLeadingZeros
isRowFromTriMatrix :: T.Matrix -> Int -> Bool
isRowFromTriMatrix mat n = hasLeadingZeros n (M.row mat n)

-- | Helper function to check if the given matrix is upper triangular
upperTriHelper :: T.Matrix -> Int -> Bool
upperTriHelper mat n = P.foldl (&&) True $ P.map (isRowFromTriMatrix mat) [0..(n-1)]

-- | to check if the given matrix is upper triangular
isUpperTri :: T.Matrix -> Bool
isUpperTri mat = upperTriHelper mat (fst (M.dimension mat))

-- | to check if the given matrix is invertible
isInvertible :: T.Matrix -> Bool
isInvertible mat = isFloatZero (det mat) 0.0000001

------------------------------- matrix manipulation functions -------------------------

-- | make all values of the matrix absolute
makeAbs :: T.Matrix -> T.Matrix
makeAbs mat = M.map P.abs mat

-- | make all values of the (row-list represented matrix)
-- | decimal values with a given precision
makeDecimal :: Word8 -> T.Matrix -> [[Decimal]]
makeDecimal prec mat =
    let
        matList = M.toRowLists mat
        rowToDecimal prec ls = P.map (realFracToDecimal prec) ls
    in
        P.map (rowToDecimal prec) matList

-- | to get the determinant of a given matrix
-- | taken from here:
-- | http://michaeljgilliland.blogspot.com/2013/04/determinant-in-haskell.html
det :: T.Matrix -> Double
det mat =
    let
        detHelper [[x]] = x
        detHelper mat = P.sum [(-1)^i*x*(detHelper (getRest i mat)) | (i, x) <- P.zip [0..] (P.head mat)]
        listMatrix = M.toRowLists mat
    in
        detHelper listMatrix

-- |  getRest function simply returns the matrix without the head row (topmost) and without the ith column
getRest :: Int -> [[Double]] -> [[Double]]
getRest col mat =
    let
        decapitatedMat = P.tail mat
    in
        P.map (leaveElement col) decapitatedMat

----------------------------- auxiliary helper functions -----------------------------------

-- | function to delete i'th element from list
leaveElement :: Ord a => Int -> [a] -> [a]
leaveElement _ [] = []
leaveElement i (x:xs)
        | i == 0    = xs
        | otherwise = x : leaveElement (i-1) xs

-- | compare a double value to some epsilon
isFloatZero :: Double -> Double -> Bool
isFloatZero n eps =
    if (abs n) > (abs eps)
        then True
    else False

max_stars_repo_name: Magalame/dh-core

max_issues_repo_name: Magalame/dh-core

max_forks_repo_name: Magalame/dh-core


text:
module LispVal
  ( LispVal(..)
  ) where

import Data.Complex

data LispVal
  = Atom String
  | List [LispVal]
  | DottedList [LispVal] LispVal
  | Number Integer
  | String String
  | Bool Bool
  | Character Char
  | Float Double
  | Ratio Rational
  | Complex (Complex Double)

max_stars_repo_name: anushriadhia/haskellScheme

max_issues_repo_name: anushriadhia/haskellScheme

max_forks_repo_name: anushriadhia/haskellScheme


text:
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import NaiveBayes as NB
import StaticTypeNet
import Mnist (readMnistAndShow, readData, dat, label, readMnist, trainData, testData, LabeledData, dataLabel, render)
import Data.Maybe
import System.Environment (getArgs)
import Control.Monad.Random (evalRandIO)
import Text.Read (readMaybe)
import Numeric.LinearAlgebra (Vector, R, toRows, fromList, toList)

main :: IO ()
main = do
  -- readMnistAndShow
  -- now read MNIST and print out the prior distribution
  mnist <- readMnist
  let trainLabeledData :: LabeledData = trainData mnist  -- train set
      testLabeledData :: LabeledData = testData mnist  -- test set
      pri = NB.prior trainLabeledData
      lk = NB.likelihood trainLabeledData
      naiveBayesClassifier :: NB.NaiveBayesClassifier = NB.NBC pri lk  -- create a classifier
      predictor = predict naiveBayesClassifier

  -- show the prior
  putStrLn $ "Below is the prior distribution"
  putStrLn $ show pri  -- show prior values

  let datlab :: [(Vector R, R)] = dataLabel testLabeledData  -- data-label pairs
      predictedAndTruth ld = (predictor (fst ld), floor $ snd ld)
      -- the result of prediction
      predictedResult :: [(Int, Int)] = map predictedAndTruth datlab
      -- calculate the accuracy
      accuracy :: [(Int, Int)] -> Double
      accuracy res = NB.fractionDiv
          (foldl
              (\acc elem -> if fst elem == snd elem
                            then acc + 1
                            else acc)
              0 res)  -- accumulate correct predictions
          (length res)  -- ...divided by total size

  -- show predicted value and truth label
  mapM_ (putStrLn . show . predictedAndTruth) datlab

  putStrLn $ "The accuracy is: "
  putStrLn $ show $ accuracy predictedResult

  -- show the 100th sample
  putStrLn $ render $ fst $ datlab !! 100
  putStrLn $ show $ predictedResult !! 100

  args <- getArgs
  let n = readMaybe =<< (args !!? 0)
      rate = readMaybe =<< (args !!? 1)
  putStrLn "Training network .."
  putStrLn =<< evalRandIO (netTest (fromMaybe 0.25 rate) (fromMaybe 500000 n))

max_stars_repo_name: deNsuh/hasml

max_issues_repo_name: deNsuh/hasml

max_forks_repo_name: deNsuh/hasml