text:
{-# LANGUAGE FlexibleContexts, ScopedTypeVariables, MultiParamTypeClasses, FlexibleInstances #-}

{-| Various utiliy signal processing functions -}
module SDR.Util (
    -- * Classes
    Mult,
    mult,

    -- * Conversion to floating point for reception
    -- ** RTLSDR
    interleavedIQUnsigned256ToFloat,
    interleavedIQUnsignedByteToFloat,
    interleavedIQUnsignedByteToFloatSSE,
    interleavedIQUnsignedByteToFloatAVX,
    interleavedIQUnsignedByteToFloatFast,

    -- ** BladeRF
    interleavedIQSigned2048ToFloat,
    interleavedIQSignedWordToFloat,
    interleavedIQSignedWordToFloatSSE,
    interleavedIQSignedWordToFloatAVX,
    interleavedIQSignedWordToFloatFast,

    -- * Conversion from floating point for transmission
    -- ** BladeRF
    complexFloatToInterleavedIQSigned2048,
    complexFloatToInterleavedIQSignedWord,

    -- * Scaling
    scaleC,
    scaleCSSE,
    scaleCAVX,
    scaleFast,

    -- * Mapping over complex numbers
    cplxMap,

    -- * Frequency shifting
    halfBandUp,
    quarterBandUp,

    -- * Data streams
    streamString,
    streamRandom,

    -- * Automatic gain control
    agc,
    agcPipe,

    -- * Squashing initialization into the Pipe
    combineInit,
    combineInitTrans
    ) where

import           Foreign.C.Types
import           Data.Complex
import qualified Data.Vector.Generic          as VG
import qualified Data.Vector.Generic.Mutable  as VGM
import qualified Data.Vector.Storable         as VS
import qualified Data.Vector.Storable.Mutable as VSM
import           Control.Monad.Primitive
import           Data.Coerce
import           Foreign.Ptr
import           System.IO.Unsafe
import           Foreign.Storable.Complex
import           Control.Monad
import qualified System.Random.MWC as R
import           Data.Bits
import           Pipes
import qualified Pipes.Prelude as P
import           Data.Word
import           Foreign.Storable
import           Control.Arrow as A
import           Data.Tuple

import           SDR.CPUID
import           SDR.VectorUtils
import           SDR.PipeUtils

-- | A class for things that can be multiplied by a scalar.
class Mult a b where
    mult :: a -> b -> a

instance (Num a) => Mult a a where
    mult = (*)

instance (Num a) => Mult (Complex a) a where
    mult (x :+ y) z = (x * z) :+ (y * z)

-- | Create a vector of complex floating samples from a vector of interleaved I Q components. Each input element ranges from 0 to 255. This is the format that RTLSDR devices use.
{-# INLINE interleavedIQUnsigned256ToFloat #-}
interleavedIQUnsigned256ToFloat :: (Num a, Integral a, Num b, Fractional b, VG.Vector v1 a, VG.Vector v2 (Complex b)) => v1 a -> v2 (Complex b)
interleavedIQUnsigned256ToFloat input = VG.generate (VG.length input `quot` 2) convert
    where
    {-# INLINE convert #-}
    convert idx  = convert' (input `VG.unsafeIndex` (2 * idx)) :+ convert' (input `VG.unsafeIndex` (2 * idx + 1))
    {-# INLINE convert' #-}
    convert' val = (fromIntegral val - 128) / 128

foreign import ccall unsafe "convertC"
    convertC_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C and specialized for unsigned byte inputs and Float outputs.
interleavedIQUnsignedByteToFloat :: VS.Vector CUChar -> VS.Vector (Complex Float)
interleavedIQUnsignedByteToFloat inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertC_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

foreign import ccall unsafe "convertCSSE"
    convertCSSE_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using SSE intrinsics and specialized for unsigned byte inputs and Float outputs.
interleavedIQUnsignedByteToFloatSSE :: VS.Vector CUChar -> VS.Vector (Complex Float)
interleavedIQUnsignedByteToFloatSSE inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertCSSE_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

foreign import ccall unsafe "convertCAVX"
    convertCAVX_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using AVX intrinsics and specialized for unsigned byte inputs and Float outputs.
interleavedIQUnsignedByteToFloatAVX :: VS.Vector CUChar -> VS.Vector (Complex Float)
interleavedIQUnsignedByteToFloatAVX inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertCAVX_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

-- | Same as `interleavedIQUnsigned256ToFloat` but uses the fastest SIMD instruction set your processor supports and specialized for unsigned byte inputs and Float outputs.
interleavedIQUnsignedByteToFloatFast :: CPUInfo -> VS.Vector CUChar -> VS.Vector (Complex Float)
interleavedIQUnsignedByteToFloatFast info = featureSelect info interleavedIQUnsignedByteToFloat [(hasAVX2, interleavedIQUnsignedByteToFloatAVX), (hasSSE42, interleavedIQUnsignedByteToFloatSSE)]

-- | Create a vector of complex float samples from a vector of interleaved I Q components. Each input element ranges from -2048 to 2047. This is the format that the BladeRF uses.
{-# INLINE interleavedIQSigned2048ToFloat #-}
interleavedIQSigned2048ToFloat :: (Num a, Integral a, Num b, Fractional b, VG.Vector v1 a, VG.Vector v2 (Complex b)) => v1 a -> v2 (Complex b)
interleavedIQSigned2048ToFloat input = VG.generate (VG.length input `quot` 2) convert
    where
    {-# INLINE convert #-}
    convert idx  = convert' (input `VG.unsafeIndex` (2 * idx)) :+ convert' (input `VG.unsafeIndex` (2 * idx + 1))
    {-# INLINE convert' #-}
    convert' val = fromIntegral val / 2048

foreign import ccall unsafe "convertCBladeRF"
    convertCBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C and specialized for signed short inputs and Float outputs.
interleavedIQSignedWordToFloat :: VS.Vector CShort -> VS.Vector (Complex Float)
interleavedIQSignedWordToFloat inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertCBladeRF_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

foreign import ccall unsafe "convertCSSEBladeRF"
    convertCSSEBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using SSE intrinsics and specialized for signed short inputs and Float outputs.
interleavedIQSignedWordToFloatSSE :: VS.Vector CShort -> VS.Vector (Complex Float)
interleavedIQSignedWordToFloatSSE inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertCSSEBladeRF_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

foreign import ccall unsafe "convertCAVXBladeRF"
    convertCAVXBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO ()

-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using AVX intrinsics and specialized for signed short inputs and Float outputs.
interleavedIQSignedWordToFloatAVX :: VS.Vector CShort -> VS.Vector (Complex Float)
interleavedIQSignedWordToFloatAVX inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf `quot` 2
    VS.unsafeWith inBuf $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            convertCAVXBladeRF_c (fromIntegral $ VG.length inBuf) iPtr oPtr
    VG.freeze outBuf

-- | Same as `interleavedIQSigned2048ToFloat` but uses the fastest SIMD instruction set your processor supports and specialized for signed short inputs and Float outputs.
interleavedIQSignedWordToFloatFast :: CPUInfo -> VS.Vector CShort -> VS.Vector (Complex Float)
interleavedIQSignedWordToFloatFast info = featureSelect info interleavedIQSignedWordToFloat [(hasAVX2, interleavedIQSignedWordToFloatAVX), (hasSSE42, interleavedIQSignedWordToFloatSSE)]

-- | Create a vector of interleaved I Q component integral samples from a vector of complex Floats. Each input ranges from -2048 to 2047. This is the format the BladeRF uses.
complexFloatToInterleavedIQSigned2048 :: (Integral b, RealFrac a, VG.Vector v1 (Complex a), VG.Vector v2 b) => v1 (Complex a) -> v2 b
complexFloatToInterleavedIQSigned2048 input = VG.generate (VG.length input * 2) convert
    where
    {-# INLINE convert #-}
    convert idx
        | even idx = convert' $ realPart (input `VG.unsafeIndex` (idx `quot` 2))
        | odd  idx = convert' $ imagPart (input `VG.unsafeIndex` (idx `quot` 2))
    {-# INLINE convert' #-}
    convert' val = round $ val * 2048

foreign import ccall unsafe "convertBladeRFTransmit"
    convertBladeRFTransmit_c :: CInt -> Ptr CFloat -> Ptr CShort -> IO ()

-- | Same as `complexFloatToInterleavedIQSigned2048` but written in C and specialized for Float inputs and signed short outputs.
complexFloatToInterleavedIQSignedWord :: VS.Vector (Complex Float) -> VS.Vector CShort
complexFloatToInterleavedIQSignedWord inBuf = unsafePerformIO $ do
    outBuf <- VGM.new $ VG.length inBuf * 2
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith outBuf $ \oPtr ->
            convertBladeRFTransmit_c (fromIntegral $ VG.length inBuf * 2) iPtr oPtr
    VG.freeze outBuf

-- | Scaling
foreign import ccall unsafe "scale"
    scale_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

-- | Scale a vector, written in C
scaleC :: Float                      -- ^ Scale factor
       -> VS.Vector Float            -- ^ Input vector
       -> VS.MVector RealWorld Float -- ^ Output vector
       -> IO ()
scaleC factor inBuf outBuf =
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            scale_c (fromIntegral (VG.length inBuf)) (coerce factor) iPtr oPtr

foreign import ccall unsafe "scaleSSE"
    scaleSSE_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat-> IO ()

-- | Scale a vector, written in C using SSE intrinsics
scaleCSSE :: Float                      -- ^ Scale factor
          -> VS.Vector Float            -- ^ Input vector
          -> VS.MVector RealWorld Float -- ^ Output vector
          -> IO ()
scaleCSSE factor inBuf outBuf =
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            scaleSSE_c (fromIntegral (VG.length inBuf)) (coerce factor) iPtr oPtr

foreign import ccall unsafe "scaleAVX"
    scaleAVX_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

-- | Scale a vector, written in C using AVX intrinsics
scaleCAVX :: Float                      -- ^ Scale factor
          -> VS.Vector Float            -- ^ Input vector
          -> VS.MVector RealWorld Float -- ^ Output vector
          -> IO ()
scaleCAVX factor inBuf outBuf =
    VS.unsafeWith (coerce inBuf) $ \iPtr ->
        VSM.unsafeWith (coerce outBuf) $ \oPtr ->
            scaleAVX_c (fromIntegral (VG.length inBuf)) (coerce factor) iPtr oPtr

-- | Scale a vector. Uses the fastest SIMD instruction set your processor supports.
scaleFast :: CPUInfo -> Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()
scaleFast info = featureSelect info scaleC [(hasAVX, scaleCAVX), (hasSSE42, scaleCSSE)]

-- | Apply a function to both parts of a complex number
cplxMap :: (a -> b)  -- ^ The function
        -> Complex a -- ^ Input complex number
        -> Complex b -- ^ Output complex number
cplxMap f (x :+ y) = f x :+ f y

-- | Multiplication by this vector shifts all frequencies up by 1/2 of the sampling frequency
halfBandUp :: (VG.Vector v n, Num n)
           => Int -- ^ The length of the Vector
           -> v n
halfBandUp size = VG.generate size func
    where
    func idx
        | even idx  = 1
        | otherwise = -1

-- | Multiplication by this vector shifts all frequencies up by 1/4 of the sampling frequency
quarterBandUp :: (VG.Vector v (Complex n), Num n)
              => Int -- ^ The length of the Vector
              -> v (Complex n)
quarterBandUp size = VG.generate size func
    where
    func idx
        | m == 0 = 1    :+ 0
        | m == 1 = 0    :+ 1
        | m == 2 = (-1) :+ 0
        | m == 3 = 0    :+ (-1)
        where
        m = idx `mod` 4

-- | A Producer that streams vectors of the bits that make up the string argument concatenated repeatedly. Each bit is encoded as a float with value (+1) for 1 and (-1) for 0.
streamString :: forall m b. (FiniteBits b, Monad m)
             => [b] -- ^ The string whose bits are to be streamed
             -> Int -- ^ The size of each streamed vector
             -> Producer (VS.Vector Float) m ()
streamString str size = P.unfoldr (return . Right . func) (str, 0)
    where
    bitsPerChar = finiteBitSize (undefined :: b)
    toFloat :: Bool -> Float
    toFloat x   = if x then 1 else (-1)
    func = vUnfoldr size funcy
        where
        funcy ([], offsetChar) = funcy (str, 0)
        funcy (rem@(x:xs), offsetChar)
            | offsetChar == bitsPerChar = funcy (xs, 0)
            | otherwise                 = (toFloat $ testBit x offsetChar, (rem, offsetChar + 1))

-- | A Producer that streams vectors of random bits. Each bit is encoded as a float with value (+1) for 1 and (-1) for 0.
streamRandom :: forall m. PrimMonad m
             => Int -- ^ The size of each streamed vector
             -> Producer (VS.Vector Float) m ()
streamRandom size = do
    gen   <- lift R.create
    start <- lift $ R.uniform gen
    P.unfoldr (liftM Right . func gen) (start, 0)
    where
    toFloat :: Bool -> Float
    toFloat x   = if x then 1 else (-1)
    func :: R.Gen (PrimState m) -> (Word64, Int) -> m (VS.Vector Float, (Word64, Int))
    func gen = vUnfoldrM size funcy
        where
        funcy (current, offset) = do
            let res =  toFloat $ testBit current offset
            if offset == 63 then do
                current' <- R.uniform gen
                return (res, (current', 0))
            else return (res, (current, offset+1))

(a :+ b) `cdiv` y = (a/y) :+ (b/y)
(a :+ b) `cmul` y = (a*y) :+ (b*y)

-- | Simple automatic gain control
agc :: (Num a, Storable a, RealFloat a)
    => a                          -- ^ a
    -> a                          -- ^ reference
    -> a                          -- ^ initial state
    -> VS.Vector (Complex a)      -- ^ input vector
    -> (a, VS.Vector (Complex a)) -- ^ (final state, output vector)
agc mu reference state input = A.first snd $ swap $ vUnfoldr (VS.length input) go (0, state)
    where
    go (offset, state) =
        let
            corrected = (input VS.! offset) `cmul` state
            state'    = state + mu * (reference - magnitude corrected)
        in  (corrected, (offset + 1, state'))

-- | Simple automatic gain control pipe
agcPipe :: (Num a, Storable a, RealFloat a, Monad m)
        => a -- ^ a
        -> a -- ^ reference
        -> Pipe (VS.Vector (Complex a)) (VS.Vector (Complex a)) m ()
agcPipe mu reference = pMapAccum (agc mu reference) 1

-- | Specializes to combineInit :: IO (Pipe a b IO ()) -> Pipe a b IO ()
combineInit :: (Monad m, MonadTrans t, Monad (t m)) => m (t m a) -> t m a
combineInit = join . lift

-- | Specializes to combineInitTrans :: EitherT String IO (Pipe a b IO ()) -> Pipe a b (EitherT String IO) ()
combineInitTrans :: (Monad (t1 m), Monad (t (t1 m)), MonadTrans t, Monad m, MFunctor t, MonadTrans t1) => (t1 m) ((t m) a) -> t (t1 m) a
combineInitTrans = combineInit . fmap (hoist lift)

max_stars_repo_name: adamwalker/sdr

max_issues_repo_name: peixian/sdr

max_forks_repo_name: adamwalker/sdr


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DeriveDataTypeable, DeriveGeneric #-}
-- |
-- Module    : Statistics.Distribution.CauchyLorentz
-- Copyright : (c) 2011 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- The Cauchy-Lorentz distribution. It's also known as Lorentz
-- distribution or Breit–Wigner distribution.
--
-- It doesn't have mean and variance.
module Statistics.Distribution.CauchyLorentz (
    CauchyDistribution
  , cauchyDistribMedian
  , cauchyDistribScale
    -- * Constructors
  , cauchyDistribution
  , cauchyDistributionE
  , standardCauchy
  ) where

import Control.Applicative
import Data.Aeson             (FromJSON(..), ToJSON, Value(..), (.:))
import Data.Binary            (Binary(..))
import Data.Maybe             (fromMaybe)
import Data.Data              (Data, Typeable)
import GHC.Generics           (Generic)
import qualified Statistics.Distribution as D
import Statistics.Internal

-- | Cauchy-Lorentz distribution.
data CauchyDistribution = CD {
    -- | Central value of Cauchy-Lorentz distribution which is its
    --   mode and median. Distribution doesn't have mean so function
    --   is named after median.
    cauchyDistribMedian :: {-# UNPACK #-} !Double
    -- | Scale parameter of Cauchy-Lorentz distribution. It's
    --   different from variance and specify half width at half
    --   maximum (HWHM).
  , cauchyDistribScale  :: {-# UNPACK #-} !Double
  }
  deriving (Eq, Typeable, Data, Generic)

instance Show CauchyDistribution where
  showsPrec i (CD m s) = defaultShow2 "cauchyDistribution" m s i
instance Read CauchyDistribution where
  readPrec = defaultReadPrecM2 "cauchyDistribution" cauchyDistributionE

instance ToJSON   CauchyDistribution
instance FromJSON CauchyDistribution where
  parseJSON (Object v) = do
    m <- v .: "cauchyDistribMedian"
    s <- v .: "cauchyDistribScale"
    maybe (fail $ errMsg m s) return $ cauchyDistributionE m s
  parseJSON _ = empty

instance Binary CauchyDistribution where
    put (CD m s) = put m >> put s
    get = do
      m <- get
      s <- get
      maybe (error $ errMsg m s) return $ cauchyDistributionE m s

-- | Cauchy distribution
cauchyDistribution :: Double    -- ^ Central point
                   -> Double    -- ^ Scale parameter (FWHM)
                   -> CauchyDistribution
cauchyDistribution m s
  = fromMaybe (error $ errMsg m s)
  $ cauchyDistributionE m s

-- | Cauchy distribution
cauchyDistributionE :: Double    -- ^ Central point
                    -> Double    -- ^ Scale parameter (FWHM)
                    -> Maybe CauchyDistribution
cauchyDistributionE m s
  | s > 0     = Just (CD m s)
  | otherwise = Nothing

errMsg :: Double -> Double -> String
errMsg _ s
  = "Statistics.Distribution.CauchyLorentz.cauchyDistribution: FWHM must be positive. Got "
  ++ show s

-- | Standard Cauchy distribution. It's centered at 0 and have 1 FWHM
standardCauchy :: CauchyDistribution
standardCauchy = CD 0 1

instance D.Distribution CauchyDistribution where
  cumulative (CD m s) x = 0.5 + atan( (x - m) / s ) / pi

instance D.ContDistr CauchyDistribution where
  density (CD m s) x = (1 / pi) / (s * (1 + y*y))
    where y = (x - m) / s
  quantile (CD m s) p
    | p > 0 && p < 1 = m + s * tan( pi * (p - 0.5) )
    | p == 0         = -1 / 0
    | p == 1         =  1 / 0
    | otherwise      =
      error $ "Statistics.Distribution.CauchyLorentz.quantile: p must be in [0,1] range. Got: "++show p
  complQuantile (CD m s) p
    | p > 0 && p < 1 = m + s * tan( pi * (0.5 - p) )
    | p == 0         =  1 / 0
    | p == 1         = -1 / 0
    | otherwise      =
      error $ "Statistics.Distribution.CauchyLorentz.complQuantile: p must be in [0,1] range. Got: "++show p

instance D.ContGen CauchyDistribution where
  genContVar = D.genContinuous

instance D.Entropy CauchyDistribution where
  entropy (CD _ s) = log s + log (4*pi)

instance D.MaybeEntropy CauchyDistribution where
  maybeEntropy = Just . D.entropy

ext: hs


text:
module Fungine.Graphics.Model where

import Fungine.Prelude
import Fungine.Graphics.Geometry
import Fungine.Graphics.Texture
import Numeric.LinearAlgebra

type BoneId = Word16
data Bone = Bone Int [Bone]

data MeshVertex = MeshVertex { vertexBones :: [(BoneId,Float)]
                             , vertex :: Point3d
                             , vertexUVCoord :: Point2d
                             , vertexNormal :: Point3d
                             }

data Model2d = Model2d { m2dRectangle :: Rectangle
                     , m2dTexture :: Texture
                     , m2dUVCoord :: Rectangle
}

data Model3d = Model3d   { m3dBones :: Map BoneId Bone
                       , m3dMesh :: Vector MeshVertex
                       , m3dTexture :: Texture
                       }

max_stars_repo_name: jmdurr/fungine

max_issues_repo_name: jmdurr/fungine

max_forks_repo_name: jmdurr/fungine


text:
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}

{-# LANGUAGE ViewPatterns #-}

import           Numeric.Sundials.ARKode.ODE
import           Numeric.LinearAlgebra
import           Graphics.Plot

vanderpol mu = do
    let xdot nu _t [x,v] = [v, -x + nu * v * (1-x*x)]
        xdot _ _ _ = error "vanderpol RHS not defined"
        ts = linspace 1000 (0,50)
        sol = toColumns $ odeSolve (xdot mu) [1,0] ts
    mplot (ts : sol)
    mplot sol

harmonic w d = do
    let xdot u dd _t [x,v] = [v, a*x + b*v] where a = -u*u; b = -2*dd*u
        xdot _ _ _ _ = error "harmonic RHS not defined"
        ts = linspace 100 (0,20)
        sol = odeSolve (xdot w d) [1,0] ts
    mplot (ts : toColumns sol)

kepler v a = mplot (take 2 $ toColumns sol) where
    xdot _t [x,y,vx,vy] = [vx,vy,x*k,y*k]
        where g=1
              k=(-g)*(x*x+y*y)**(-1.5)
    xdot _ _ = error "kepler RHS not defined"
    ts = linspace 100 (0,30)
    sol = odeSolve xdot [4, 0, v * cos (a*degree), v * sin (a*degree)] ts
    degree = pi/180

main = do
    vanderpol 2
    harmonic 1 0
    harmonic 1 0.1
    kepler 0.3 60
    kepler 0.4 70
    vanderpol' 2

-- example of odeSolveV with jacobian
vanderpol' mu = do
    let xdot nu _t (toList->[x,v]) = fromList [v, -x + nu * v * (1-x*x)]
        xdot _ _ _ = error "vanderpol' RHS not defined"
        jac _ (toList->[x,v]) = (2><2) [ 0          ,          1
                                       , -1-2*x*v*mu, mu*(1-x**2) ]
        jac _ _ = error "vanderpol' Jacobian not defined"
        ts = linspace 1000 (0,50)
        hi = pure $ (ts!1 - ts!0) / 100.0
        sol = toColumns $ odeSolveV (SDIRK_5_3_4 jac) hi 1E-8 1E-8 (xdot mu) (fromList [1,0]) ts
    mplot sol

max_stars_repo_name: regnat/hmatrix-sundials

max_issues_repo_name: regnat/hmatrix-sundials

max_forks_repo_name: regnat/hmatrix-sundials


text:
module NeuralNet.Train (
  nnBackward,
  updateNNParams
) where

import NeuralNet.Net
import NeuralNet.Layer
import NeuralNet.Example
import NeuralNet.Activation
import NeuralNet.Matrix
import Numeric.LinearAlgebra

nnBackward :: NeuralNet -> [ForwardPropStep] -> ExampleSet -> [(Matrix Double, Matrix Double)]
nnBackward nn steps examples
  | numLayers /= stepsSize - 1 = error ("should be one entry (" ++ show stepsSize ++ ") for each layer (" ++ show numLayers ++ ")")
  | otherwise                  = reverse (nnBackwardStep layers orderedSteps dAL)
    where
      layers = reverse (nnLayers nn)
      orderedSteps = reverse steps
      numLayers = length layers
      stepsSize = length steps

      y = exampleSetY examples
      al = forwardPropA (head orderedSteps)
      unitRowVector = row (replicate (cols y) 1)
      dAL1 = y / al
      dAL2 = (unitRowVector - y) / (unitRowVector - al)
      dAL = -(dAL1 - dAL2)

nnBackwardStep :: [NeuronLayer] -> [ForwardPropStep] -> Matrix Double -> [(Matrix Double, Matrix Double)]
nnBackwardStep [] _ _ = []
nnBackwardStep _ [] _ = error "Shouldn't get here"
nnBackwardStep _ [_] _ = error "Shouldn't get here"
nnBackwardStep (l:ls) (s:ss@(ps:_)) prevDA = (dW, db) : nnBackwardStep ls ss dA
  where
    (dA, dW, db) = nnBackwardLayer l prevDA (forwardPropA ps) (forwardPropZ s)

nnBackwardLayer :: NeuronLayer -> Matrix Double -> Matrix Double -> Matrix Double -> (Matrix Double, Matrix Double, Matrix Double)
nnBackwardLayer l dA prevA z = (dANext, dW, db)
  where
    m = fromIntegral (cols prevA)
    dZ = backward (layerActivation l) dA z
    dW = scale (1/m) (dZ <> tr prevA)
    db = scale (1/m) (sumRows dZ)
    dANext = tr (layerW l) <> dZ

updateNNParams :: NeuralNet -> [(Matrix Double, Matrix Double)] -> Double -> NeuralNet
updateNNParams nn grads learningRate
  | length layers /= length grads = error "Must provide one grad per layer"
  | otherwise                     = updateNNLayers nn newLayers
    where
      layers = nnLayers nn
      layerSets = zip layers grads
      newLayers = map (\(l, (dW, dB)) -> updateLayerParamsFromGrads l dW dB learningRate) layerSets

updateLayerParamsFromGrads :: NeuronLayer -> Matrix Double -> Matrix Double -> Double -> NeuronLayer
updateLayerParamsFromGrads l dW db learningRate = updateLayerParams l newDw newDb
  where
    newDw = layerW l - scale learningRate dW
    newDb = layerB l - scale learningRate db

max_stars_repo_name: danielholmes/neural-net

max_issues_repo_name: danielholmes/neural-net

max_forks_repo_name: danielholmes/neural-net


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UnicodeSyntax #-}
module Day20 ( runSolution) where

import Data.ByteString (ByteString)
import qualified Data.ByteString as BS
import Data.Attoparsec.ByteString (Parser)
import qualified Data.Attoparsec.ByteString as DAB
import Data.Attoparsec.ByteString.Char8 (anyChar, space, notChar, decimal, string, char, endOfLine)

import Data.Array.IArray (Array, (!))
import qualified Data.Array.IArray as IA
import qualified Data.Complex as DC
import qualified Linear as L
import Control.Arrow
import qualified Data.List as DL
import Control.Lens
import Data.Function (on)
import Data.Maybe (fromJust)
import qualified Control.Monad as CM
import qualified Data.Set as Set
import qualified Text.PrettyPrint.Boxes as Boxes
import qualified Control.Foldl as L

import Debug.Trace

type Picture = Array (L.V2 Int) Char

type Piece = (Int, Picture)

-- A `FittingRow` is a list of pieces that fit side by side in an horizontal
-- line.
type FittingRow = [Piece]

-- A `FittingPuzzle` is a list of rows where each row matches the one above and
-- below.
type FittingPuzzle = [FittingRow]

data HorizontalGrowthDirection = ToTheLeft | ToTheRight deriving Show
data VerticalGrowthDirection = Upwards | Downwards deriving Show

-- | Build a solution adding pieces to the right
--
-- For example, starting with P0 it adds P1 onwards where
-- P1 matches the *right* edge of P0.
--
-- (P0) P1 P2 P3
--
-- It returns all posibilities. That is
--
-- [] : Don't add more pieces
-- [P1] : Just add P1
-- [P1, P2] : Just add P1 and P2
-- ...
--
-- This leaves free pieces to build other parts of the solution. E.g: to the
-- left of P0.
--
-- It doesn't include the starting piece (P0)
--
-- if HorizontalGrowthDirection is ToTheLeft. P1 fits to the left of P0
--
-- P1 (P0)
--
-- regardless the resulting rows are in reverse order. That is for:
--
-- P3 P2 P1 (P0)
--
-- the output is [[],[P1],[P1,P2],[P1,P2,P3]]
--
-- This means that in this case the alternatives shouls be reversed to be
-- displayed in the correct order.
topHorizontalSolution :: HorizontalGrowthDirection -> Piece -> [(Int, [Picture])] -> [FittingRow]
topHorizontalSolution growthDir (_, topLeftPic) pieces = unfoldrM unfolder (topLeftPic, pieces)
	where
	unfolder :: (Picture, [(Int, [Picture])]) -> [Maybe (Piece, (Picture, [(Int, [Picture])]))]
	-- `refPic` is the picture the next piece should fit on the side.
	-- `pieces` are the remmaining unused pieces
	--
	-- The alternatives are:
	-- 	Stop
	-- 	One alternative for each piece that matches the refPic on the side.
	unfolder (refPic, pieces) = (Nothing :) $
		fmap
			(\sidePiece -> Just (sidePiece,
				( snd sidePiece -- The next refPic
				, removePiece sidePiece pieces))) -- The next pieces
			(concatMap filterSideMatchs pieces)
		where

		filterSideMatchs :: (Int, [Picture]) -> [(Int, Picture)]
		filterSideMatchs (pid, trans) = (pid,) <$> filter selector trans

		selector = case growthDir of
			ToTheRight -> match Sideways refPic
			ToTheLeft -> (\t -> match Sideways t refPic)

-- | Build a solution adding pieces to both sides
topMiddleSolution :: (Int, Picture) -> [(Int, [Picture])] -> [FittingRow]
topMiddleSolution topCenter@(_, topCenterPic) pieces = [topCenter] : do
	rightPiece <- pieces >>= filterRightMatches
	leftPiece <- removePiece rightPiece pieces >>= filterLeftMatches

	let remPieces = removePiece rightPiece (removePiece leftPiece pieces)

	rightSolution <- topHorizontalSolution ToTheRight rightPiece remPieces

	-- The soultions need to be reversed because they are to the left of the
	-- `topCenter` piece.
	leftSolution <- fmap reverse $ topHorizontalSolution ToTheLeft leftPiece $
		DL.foldl' (flip removePiece) remPieces rightSolution

	pure (leftSolution ++ [leftPiece, topCenter, rightPiece] ++ rightSolution)
	where
	filterRightMatches :: (Int, [Picture]) -> [(Int, Picture)]
	filterRightMatches (pid, trans) = (pid,) <$> filter (match Sideways topCenterPic) trans

	filterLeftMatches :: (Int, [Picture]) -> [(Int, Picture)]
	filterLeftMatches (pid, trans) = (pid,) <$> filter (\t -> match Sideways t topCenterPic) trans

-- | Build a solution adding a new row at the bottom
downwardSolution :: FittingRow -> [(Int, [Picture])] -> [FittingPuzzle]
downwardSolution _ [] = [[]]
downwardSolution (t : ts) pieces@(_ : _) = [] : do
	bottomPiece <- pieces >>= filterBelowMatches
	bsol <- rowSolution Downwards bottomPiece ts
		(filter (not . null . snd) $ removePiece bottomPiece pieces)
	CM.guard (length (bottomPiece : bsol) == 12)
	[bottomPiece : bsol] : do
		let remPieces = DL.foldl' (flip removePiece) pieces (bottomPiece : bsol)
		bsol' <- downwardSolution (bottomPiece : bsol) remPieces
		[ (bottomPiece : bsol) : bsol']
	where

	filterBelowMatches (pid, trans) = (pid,) <$> filter (match Below (snd t)) trans

-- | Build a row at the bottom or top of another.
--
-- Given a row of pieces (that fit together, although this funtion doesn't check
-- that this is the case)
--
-- T1 T2 T3 T4
--
-- And a piece P0 (that could match a possible T0 to the left of T1). It builds
-- all possible solutions (S1..) that fit to the right and below the T's.
--
--    T1 T2 T3 T4
-- P0 S1 S2 S3 S4
--
-- So S1 fits P0 from the right and T1 from below. Similarly T2, etc...
--
-- This function is not greedy. It returns all possible solutions. In the
-- running example that is:
--
-- [[], [S1], [S1, S2], [S1, S2, S3], [S1, S2, S3, S4]]
--
-- If there are more solutions, thouse too should be included on the list.
--
-- The solution is built to the right.
--
-- In case the VerticalGrowthDirection is Upwards the above example becomes:
--
-- P0 S1 S2 S3 S4
--    T1 T2 T3 T4
--
-- This function can be enhanced to build solutions to the left.
rowSolution :: VerticalGrowthDirection -> Piece -> FittingRow -> [(Int, [Picture])] -> [FittingRow]
rowSolution growthDir (_, sidePic) topRow pieces = unfoldrM unfolder (sidePic, (topRow, pieces))
	where
	unfolder :: (Picture, (FittingRow, [(Int, [Picture])])) ->
		[Maybe (Piece, (Picture, (FittingRow, [(Int, [Picture])])))]
	unfolder (_, ([], _)) = [Nothing]
	-- `refPic` is P0 in the comment of the parent function. It changes to (S1,
	-- S2,...) as the computation advances.
	--
	-- `(t : ts)` is the top (T1, T2, ...). It is consumed one by one.
	--
	-- `pieces` are the available pieces. The ones used are removed from the
	-- list for the next iteration.
	unfolder (refPic, (t : ts, pieces)) = (Nothing :) $
		fmap
			(\sidePiece -> Just (sidePiece,
				( snd sidePiece
				,	( ts
					, removePiece sidePiece pieces))))
			(concatMap filterBelowMatches pieces)
		where

		filterBelowMatches (pid, trans) = (pid,) <$> filter selector trans

		-- True if `p` matches P0 from the right and T1 from below/above.
		selector :: Picture -> Bool
		selector p = match Sideways refPic p && (
			case growthDir of
				Downwards -> match Below (snd t) p
				Upwards -> match Below p (snd t) )

-- | Build a solution adding a new row on the top
upwardSolution :: FittingRow -> [(Int, [Picture])] -> [FittingPuzzle]
upwardSolution _ [] = [[]]
upwardSolution (b : bs) pieces@(_ : _) = [] : do
	topPiece <- pieces >>= filterUpwardsMatches
	tsol <- rowSolution Upwards topPiece bs
		(filter (not . null . snd) $ removePiece topPiece pieces)
	CM.guard (length (topPiece : tsol) == 12)
	[topPiece : tsol] : do
		let remPieces = DL.foldl' (flip removePiece) pieces (topPiece : tsol)
		tsol' <- upwardSolution (topPiece : tsol) remPieces
		[tsol' ++ [topPiece : tsol]]
	where

	filterUpwardsMatches (pid, trans) = (pid,) <$> filter (\t -> match Below t (snd b)) trans

middleRowSolution :: FittingRow -> [(Int, [Picture])] -> [FittingPuzzle]
middleRowSolution middleRow pieces = [middleRow] : do
	ts <- upwardSolution middleRow pieces
	let remPieces = DL.foldl' (flip removePiece) pieces (concat ts)
	bs <- downwardSolution middleRow remPieces
	pure $ ts ++ [middleRow] ++ bs

removePiece :: Piece -> [(Int, [Picture])] -> [(Int, [Picture])]
removePiece piece = filter ((/= fst piece) . fst)

solver input = case transformed of
	[] -> []
	(p : ps) -> do
		trans <- snd p
		topRow <- topMiddleSolution (fst p, trans) ps
		CM.guard (length topRow == 12)
		let remPieces = DL.foldl' (flip removePiece) (p : ps) topRow
		bs <- middleRowSolution topRow remPieces
		CM.guard (length bs == 12)
		pure bs

	where
	transformed = fmap (second transforms) input

data RelPos = Below | Sideways deriving Show

match :: RelPos -> Picture -> Picture -> Bool
match relPos refPic testPic = and $ case relPos of
	Below -> zipWith (==)
		(fmap ( (refPic !) . (`L.V2` rows)) [1..cols]) -- the bottom row of refPic
		(fmap ( (testPic !) . (`L.V2` 1)) [1..cols]) -- the top row of testPic
	Sideways -> zipWith (==)
		(fmap ( (refPic !) . (cols `L.V2`)) [1..rows]) -- the rightmost col of refPic
		(fmap ( (testPic !) . (1 `L.V2`)) [1..rows]) -- the leftmost col of testPic
	where
	(_, L.V2 cols rows) = IA.bounds refPic

-- All unique transformations of the picture
transforms :: Picture -> [Picture]
transforms basePic =
	-- try to exploit the picture's symmetry. Flip symmetries.
	Set.toList $ Set.fromList $
	-- An Horizontal flip is a Vertical flip and a 180 degree rotation.
	rotations ++ ( (Vertical 🤸) <$> rotations )
	where
	-- try to exploit the picture's symmetry. Stop ASAP. 90 degree rotation or
	-- 180 degree rotation symmetry
	rotations = basePic : takeWhile (/= basePic) (tail $ iterate (ClockWise 🥏) basePic )

-- | Removes the edges from a picture.
removeEdge :: Picture -> Picture
removeEdge pic = IA.array (L.V2 1 1, L.V2 (xlim - 2) (ylim - 2)) $
		zip (fmap (+ L.V2 (-1) (-1)) nonBorderPixels) (fmap (pic !) nonBorderPixels)
	where
	nonBorderPixels = L.V2 <$> [(orgx + 1)..(xlim-1)] <*> [(orgy + 1)..(ylim-1)]
	(L.V2 orgx orgy, L.V2 xlim ylim) = IA.bounds pic

-- | Construct one big picture from a fitting puzzle
bigPicture :: FittingPuzzle -> Picture
bigPicture pzl = IA.array (L.V2 1 1, L.V2 xlim ylim) $ concatMap (\(y, cs) ->
	fmap (first (`L.V2` y)) cs)  $ zip [1..] $ fmap (zip [1..]) picStrs
	where
	ylim = length picStrs
	xlim = length $ head picStrs
	-- Fold all rows together
	picStrs = concatMap rows pictures
	-- Build a big-row
	rows :: [Picture] -> [String]
	rows pcs = DL.foldl' rowFolder [] pcs
	rowFolder :: [String] -> Picture -> [String]
	rowFolder [] p = pl p
	rowFolder acc@(_ : _) p = zipWith (++) acc (pl p)
	-- pictures without borders
	pictures :: [[Picture]]
	pictures = fmap (fmap (removeEdge . snd)) pzl
	-- picture lines
	pl :: Picture -> [String]
	pl p = fmap (fmap snd) $ DL.groupBy ((==) `on` ((^. _2) . fst)) $ DL.sortOn ((^. _2) . fst) $  IA.assocs p

-- | Parse the input into a 2D map of chars using a default character if some
-- lines are shorter
parsePicture :: Char -> Parser Picture
parsePicture def = do
	parsedLines <- parseLines
	let maxRowLength =  (maximum . fmap length) parsedLines
	let nRows = length parsedLines
	pure $ (IA.array (L.V2 1 1, L.V2 maxRowLength nRows) . concatMap mapper . zip [1..nRows] . fmap (makeRow maxRowLength)) parsedLines
	where
	mapper :: (Int, [(Int, Char)]) -> [(L.V2 Int, Char)]
	mapper (y, row) = fmap (first (`L.V2` y)) row
	parseLines :: Parser [[Char]]
	parseLines = parseLine `DAB.sepBy1'` endOfLine
	parseLine :: Parser [Char]
	parseLine = DAB.many1' (notChar '\n')
	makeRow :: Int -> [Char] -> [(Int, Char)]
	makeRow size content = zip [1..size] (content ++ repeat def)

parseInput :: Parser [Piece]
parseInput = DAB.many1' parseSingle
	where
	parseSingle = (,) <$>
		(string "Tile " *> decimal <* char ':' <* endOfLine) <*>
		parsePicture 'X' <* DAB.choice [endOfLine <* endOfLine, endOfLine]

-- | Cound the #'s present in a picture
-- Right answer is 1607
waterRoughness :: Picture -> Int
waterRoughness = length . filter ((=='#') . snd) . IA.assocs

-- | Detect the monsters and mark them with 'O'
detectMonsters :: Picture -> Picture
detectMonsters bigPic = bigPic IA.// fmap (,'O') monsters
	where
	-- The positions of the parts of the monsters
	monsters = concatMap (fmap fst) $ filter (all ((=='#') . snd))  $ fmap mapper origins
	mapper origin = fmap ( (id &&& (bigPic !)) . (+ origin)) monsterMold
	-- All valid positions of the mosnter mold in the bigPic
	origins = L.V2 <$> [picMinX..picMaxX - maxx] <*> [picMinY..picMaxY - maxy]
	(L.V2 picMinX picMinY, L.V2 picMaxX picMaxY) = IA.bounds bigPic
	-- The width and height of the moster mold
	(Just maxx, Just maxy) = L.fold ((,) <$> L.premap (^. L._x) L.maximum <*> L.premap (^. L._y) L.maximum) monsterMold
	-- The positions where the seaMonsterInput has a '#'
	-- Starts at 0, 0 becase it's intended to be uses as an offset from an
	-- origin.
	monsterMold :: [L.V2 Int]
	monsterMold = concatMap (\(y, xs) -> fmap ((`L.V2` y) . fst) xs ) $
		zip [0..] $
		fmap
			( filter ((== '#') . snd) . zip [0..] )
			seaMonsterInput

runSolution :: FilePath -> IO ()
runSolution filePath = do
	putStrLn "**Day 20**"
	contents <- BS.readFile filePath
	let parseResult = DAB.parseOnly parseInput contents
	case parseResult of
		Left err -> putStrLn $ "Error :" ++ err
		Right input -> do
			CM.forM_ (take 1 $ solver input) $ \ pzl -> do
				putStrLn $ Boxes.render $ picLines $ bigPicture pzl
				CM.forM_ (transforms $ bigPicture pzl) $ print . waterRoughness . detectMonsters
				putStrLn "-- End Solution --"
	putStrLn "Done"

seaMonsterInput :: [String]
seaMonsterInput =
 [ "                  # "
 , "#    ##    ##    ###"
 , " #  #  #  #  #  #   " ]

puzzleLines :: FittingPuzzle -> Boxes.Box
puzzleLines pzl = Boxes.vsep 0 Boxes.left pls
	where
	pls = fmap (Boxes.hsep 1 Boxes.top) $ fmap pieceLines <$> pzl

pieceLines :: Piece -> Boxes.Box
pieceLines (pid, pic) = Boxes.text (show pid) Boxes.// picLines pic

picLines :: Picture -> Boxes.Box
picLines pic = Boxes.vcat Boxes.left $ fmap (Boxes.text . (snd <$>)) $
	DL.groupBy ((==) `on` ((^. _2) . fst)) $ DL.sortOn ((^. _2) . fst) $  IA.assocs pic

data RotDir = ClockWise | CounterClockWise deriving (Show, Eq)

--`#rotp`
(🥏)  :: RotDir -> Picture -> Picture
(🥏)  rotDir pic = IA.array (org, L.V2 ylim xlim) $
	first (untranslate . rotator rotDir . translate) <$> IA.assocs pic
	where

	translate :: L.V2 Int -> L.V2 Float
	translate (L.V2 x y) = L.V2 (fromIntegral x - 1 - fxlim) (fromIntegral y - 1 - fylim)

	untranslate :: L.V2 Float -> L.V2 Int
	untranslate (L.V2 x y) = round <$> L.V2 (x + 1 + fxlim) (y + 1 + fylim)

	rotator :: RotDir -> L.V2 Float -> L.V2 Float
	rotator CounterClockWise (L.V2 x y) = L.V2 y (negate x)
	rotator ClockWise (L.V2 x y) = L.V2 (negate y) x

	(org, L.V2 xlim ylim) = IA.bounds pic
	fxlim :: Float
	fylim :: Float
	(fxlim, fylim) = (fromIntegral (xlim - 1) / 2, fromIntegral (ylim - 1) / 2)

data FlipAxis = Vertical | Horizontal deriving Show

--`#flpp`
(🤸) :: FlipAxis -> Picture -> Picture
(🤸) flipAxis pic = IA.array (org, L.V2 ylim xlim) $
	first (untranslate . flipper flipAxis . translate) <$> IA.assocs pic

	where

	translate :: L.V2 Int -> L.V2 Float
	translate (L.V2 x y) = L.V2 (fromIntegral x - 1 - fxlim) (fromIntegral y - 1 - fylim)

	untranslate :: L.V2 Float -> L.V2 Int
	untranslate (L.V2 x y) = round <$> L.V2 (x + 1 + fxlim) (y + 1 + fylim)

	flipper :: FlipAxis -> L.V2 Float -> L.V2 Float
	flipper Horizontal (L.V2 x y) = L.V2 x (negate y)
	flipper Vertical (L.V2 x y) = L.V2 (negate x) y

	(org, L.V2 xlim ylim) = IA.bounds pic
	fxlim :: Float
	fylim :: Float
	(fxlim, fylim) = (fromIntegral (xlim - 1) / 2, fromIntegral (ylim - 1) / 2)

-- | unfoldr in a monad.
--
-- This is missing from Data.List
unfoldrM :: Monad m => (a -> m (Maybe (b, a))) -> a -> m [b]
unfoldrM f a = f a >>= \case
	Nothing -> pure []
	Just (b, a') -> (b :) <$> unfoldrM f a'


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RecordWildCards  #-}

module Network.Layer
( LayerDefinition(..)
, Layer(..)
, Connectivity
, RandomTransform

, showableToLayer

, createLayer
, scaleLayer
, connectFully

, randomList
, boxMuller
, normals
, uniforms
, boundedUniforms
) where

import           Data.Binary           (Binary (..), decode, encode)
import           Network.Neuron
import           Numeric.LinearAlgebra
import           System.Random

-- | The LayerDefinition type is an intermediate type initialized by the
--   library user to define the different layers of the network.
data LayerDefinition = LayerDefinition { neuronDef   :: Neuron
                                       , neuronCount :: Int
                                       , connect     :: Connectivity
                                       }

-- | The Layer type, which stores the weight matrix, the bias matrix, and
--   a neuron type.
data Layer = Layer { weightMatrix :: Matrix Double
                   , biasVector   :: Vector Double
                   , neuron       :: Neuron
                   } deriving Show

instance Binary (Layer) where
  put Layer{..} = do put weightMatrix; put biasVector
  get = do weightMatrix <- get; biasVector <- get; return Layer{..}

-- | Connectivity is the type alias for a function that defines the connective
--   matrix for two layers (fully connected, convolutionally connected, etc.)
type Connectivity = Int -> Int -> Matrix Double

-- | A random transformation type alias. It is a transformation defined on an
--   infinite list of uniformly distributed random numbers, and returns a list
--   distributed on the transforming distribution.
type RandomTransform = [Double] -> [Double]

-- | The createLayer function takes in a random transformation on an infinite
--   stream of uniformly generated numbers, a source of entropy, and two
--   layer definitions, one for the previous layer and one for the next layer.
--   It returns a layer defined by the Layer type -- a weight matrix, a bias
--   vector, and a neuron type.
createLayer :: (RandomGen g)
  => RandomTransform -> g -> LayerDefinition -> LayerDefinition -> Layer
createLayer t g layerDef layerDef' =
  Layer (randomMatrix * (connectivity i j))
        (randomVector * bias)
        (neuronDef layerDef)
  where randomMatrix = (i >< j) (randomList t g')
        randomVector = i |> (randomList t g'')
        i = neuronCount layerDef'
        j = neuronCount layerDef
        connectivity = connect layerDef'
        bias = i |> (repeat 1) -- bias connectivity (full)
        (g', g'') = split g

scaleLayer :: Double -> Layer -> Layer
scaleLayer factor l =
  Layer (factor `scale` (weightMatrix l)) (factor `scale` (biasVector l)) (neuron l)

-- | The connectFully function takes the number of input neurons for a layer, i,
--   and the number of output neurons of a layer, j, and returns an i x j
--   connectivity matrix for a fully connected network.
connectFully :: Int -> Int -> Matrix Double
connectFully i j = (i >< j) (repeat 1)

-- | To go from a showable to a layer, we also need a neuron type,
--   which is an unfortunate restriction owed to Haskell's inability to
--   serialize functions.
showableToLayer :: (Layer, LayerDefinition) -> Layer
showableToLayer (s, d) = Layer (weightMatrix s) (biasVector s) (neuronDef d)

-- | Initialize an infinite random list given a random transform and a source
--   of entroy.
randomList :: RandomGen g => RandomTransform -> g -> [Double]
randomList transform = transform . randoms

-- | Define a transformation on the uniform distribution to generate
--   normally distributed numbers in Haskell (the Box-Muller transform)
boxMuller :: Double -> Double -> (Double, Double)
boxMuller x1 x2 = (z1, z2)
  where z1 = sqrt ((-2) * log x1) * cos (2 * pi * x2)
        z2 = sqrt ((-2) * log x1) * sin (2 * pi * x2)

-- | This is a function of type RandomTransform that transforms a list of
--   uniformly distributed numbers to a list of normally distributed numbers.
normals :: [Double] -> [Double]
normals (x1:x2:xs) = z1:z2:(normals xs)
  where (z1, z2) = boxMuller x1 x2
normals _ = []

-- | A non-transformation to return a list of uniformly distributed numbers
--   from a list of uniformly distributed numbers. It's really a matter of
--   naming consistency. It generates numbers on the range (0, 1]
uniforms :: [Double] -> [Double]
uniforms xs = xs

-- | An affine transformation to return a list of uniforms on the range
--   (a, b]
boundedUniforms :: (Double, Double) -> [Double] -> [Double]
boundedUniforms (lower, upper) xs = map affine xs
  where affine x = lower + x * (upper - lower)

max_stars_repo_name: mckeankylej/hwRecog

max_issues_repo_name: mckeankylej/hwRecog

max_forks_repo_name: mckeankylej/hwRecog


text:
module Statistics
    ( runJSONFold
    , receivedCreatedF
    , module Statistics.Block
    , module Statistics.Chart
    , module Statistics.CSV
    , module Statistics.Focus
    , module Statistics.Graph
    , module Statistics.MemPool
    , module Statistics.Report
    , module Statistics.Throughput
    , module Statistics.Tx
    ) where

import           Control.Foldl (Fold (..))
import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M

import           JSONLog
import           Statistics.Block
import           Statistics.Chart
import           Statistics.CSV
import           Statistics.Focus
import           Statistics.Graph
import           Statistics.MemPool
import           Statistics.Report
-- import Statistics.Relay
import           Statistics.Throughput
import           Statistics.Tx
import           Types
import           Universum
import           Util.Pipes (fold')

runJSONFold :: FilePath -> Fold IndexedJLTimedEvent a -> IO a
runJSONFold logDir fd = runParseLogs logDir $ fold' fd

receivedCreatedF :: Fold IndexedJLTimedEvent (Map TxHash (Maybe Timestamp))
receivedCreatedF = f <$> txFirstReceivedF <*> inBlockChainF
  where
    f :: Map TxHash Timestamp -> Map TxHash Timestamp -> Map TxHash (Maybe Timestamp)
    f rm cm = M.mapWithKey g rm
      where
        g :: TxHash -> Timestamp -> Maybe Timestamp
        g tx ts = maybe Nothing (\ts' -> Just $ ts' - ts) $ M.lookup tx cm


text:
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE TypeOperators       #-}

{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
module Test.Grenade.Layers.Internal.Pooling where

import           Grenade.Layers.Internal.Pooling

import           Numeric.LinearAlgebra hiding (uniformSample, konst, (===))

import           Hedgehog
import qualified Hedgehog.Gen as Gen
import qualified Hedgehog.Range as Range

import qualified Test.Grenade.Layers.Internal.Reference as Reference
import           Test.Hedgehog.Compat

prop_poolForwards_poolBackwards_behaves_as_reference =
  let ok extent kernel = [stride | stride <- [1..extent], (extent - kernel) `mod` stride == 0]
      output extent kernel stride = (extent - kernel) `div` stride + 1
  in  property $ do
        height   <- forAll $ choose 2 100
        width    <- forAll $ choose 2 100
        kernel_h <- forAll $ choose 1 (height - 1)
        kernel_w <- forAll $ choose 1 (width - 1)
        stride_h <- forAll $ Gen.element (ok height kernel_h)
        stride_w <- forAll $ Gen.element (ok width kernel_w)
        input    <- forAll $ (height >< width) <$> Gen.list (Range.singleton $ height * width) (Gen.realFloat $ Range.linearFracFrom 0 (-100) 100)

        let outFast       = poolForward 1 height width kernel_h kernel_w stride_h stride_w input
        let retFast       = poolBackward 1 height width kernel_h kernel_w stride_h stride_w input outFast

        let outReference  = Reference.poolForward kernel_h kernel_w stride_h stride_w (output height kernel_h stride_h) (output width kernel_w stride_w) input
        let retReference  = Reference.poolBackward kernel_h kernel_w stride_h stride_w  input outReference

        outFast === outReference
        retFast === retReference

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

max_issues_repo_name: Alien-Inc/grenade

max_forks_repo_name: Alien-Inc/grenade


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

import           Control.Applicative
import           Control.Arrow                 ((&&&))
import           Control.Lens                  hiding ((|>))
import           Control.Monad                 (liftM,void)
import           Control.Monad.Primitive.Class
import           Data.Csv                      (HasHeader(..))
import           Data.List                     (mapAccumL)
import           Data.Ratio
import           Data.Time
import qualified Data.Vector                   as V (Vector)
import qualified Data.Vector.Generic           as V hiding (Vector)
import           Graphics.EasyPlot
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Util
import           Pipes
import qualified Pipes.ByteString              as PB
import           Pipes.Csv                     (decode)
import           Pipes.Safe                    (SafeT, runSafeT)
import qualified Pipes.Safe.Prelude            as PS
import qualified Pipes.Vector                  as PV
import           System.IO                     (IOMode (ReadMode))

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

main :: IO ()
main = do
  runLRStrategy "Cumulative sum of " absReturn cmlAbsReturn 0
  runLRStrategy "Cumulative product of " pctReturn cmlPctReturn 0
  runKFStrategy "Cumulative sum of " absReturn cmlAbsReturn 0

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

-- Linear regression strategy inspired by http://businessforecastblog.com/predicting-the-sp-500-or-the-spy-exchange-traded-fund/

runLRStrategy :: String -> Return -> Return -> Matrix Double -> IO ()
runLRStrategy titleprefix freturn fcmlreturn minret = do
  (mx,my,vdtd) <- prepareData 30 freturn
  let
    t = 2484 -- number of training data points
    mm = takeRows t mx <\> takeRows t my
    my' = mx <> mm
    mr = step (my' - minret) * my
    vy = head $ toColumns my
    vr = head $ toColumns mr
    vcy = fcmlreturn vy
    vcr = fcmlreturn vr
    lcy = V.toList vcy
    lcr = V.toList vcr
    ldt = map (fromRational . dayToYear . readDay) $ V.toList vdtd
  void $ plot' [] X11 [
    Data2D [Title $ titleprefix ++ "in sample strategy returns", Style Lines] [] $ take t $ zip ldt lcr,
    Data2D [Title $ titleprefix ++ "out of sample strategy returns", Style Lines] [] $ drop t $ zip ldt lcr,
    Data2D [Title $ titleprefix ++ "base returns", Style Lines] [] $ zip ldt lcy
    ]

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

-- Kalman filter strategy

runKFStrategy :: String -> Return -> Return -> Matrix Double -> IO ()
runKFStrategy titleprefix freturn fcmlreturn minret = do
  (mx,my,vdtd) <- prepareData 30 freturn
  let
    f = diag $ 61 |> repeat 1 :: Matrix Double
    q = diag $ 61 |> repeat 0.2 :: Matrix Double
    r = (1><1) [0.1] :: Matrix Double
    x0 = 61 |> repeat 0
    p0 = q
    s0 = State x0 p0
    (s',y') = mapAccumL (\st (xt,yt) -> (kalman (System f (fromRows [xt]) q r) st yt,(xt * (sX st)) @> 1)) s0 $ zip (toRows mx) (toRows my)
    my' = fromColumns [fromList y']
    mr = (step (my' - minret)) * my
    vy = head $ toColumns my
    vr = head $ toColumns mr
    vcy = fcmlreturn vy
    vcr = fcmlreturn vr
    lcy = V.toList vcy
    lcr = V.toList vcr
    ldt = map (fromRational . dayToYear . readDay) $ V.toList vdtd
  print $ sX s'
  void $ plot' [] X11 [
    Data2D [Title $ titleprefix ++ "strategy returns", Style Lines] [] $ zip ldt lcr,
    Data2D [Title $ titleprefix ++ "base returns", Style Lines] [] $ zip ldt lcy
    ]

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

-- Kalman code from hmatrix examples

data System = System {kF, kH, kQ, kR :: Matrix Double}
data State = State {sX :: Vector Double , sP :: Matrix Double} deriving Show
type Measurement = Vector Double

kalman :: System -> State -> Measurement -> State
kalman (System f h q r) (State x p) z = State x' p' where
    px = f <> x                            -- prediction
    pq = f <> p <> trans f + q             -- its covariance
    y  = z - h <> px                       -- residue
    cy = h <> pq <> trans h + r            -- its covariance
    k  = pq <> trans h <> inv cy           -- kalman gain
    x' = px + k <> y                       -- new state
    p' = (ident (dim x) - k <> h) <> pq    -- its covariance

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

--TODO: streamline inefficient conversion between lists and vectors

prepareData :: Int -> Return -> IO (Matrix Double,Matrix Double,V.Vector String)
prepareData g freturn = do
  vsp' <- V.map ((^._1) &&& (^._7)) `liftM` readCSV "gspc.csv"
  vvx' <- V.map ((^._1) &&& (^._7)) `liftM` readCSV "vix.csv"
  let
    vm' = V.fromList $ mergeEqual (V.toList vsp') (V.toList vvx') :: V.Vector (String,Double,Double)
    vdtd = V.map (^._1) (V.drop g vm') :: V.Vector String
    [vsp'',vvx''] = [V.map] <*> [(^._2),(^._3)] <*> [vm']
    [vsp,vvx] = map V.convert [vsp'',vvx''] :: [Vector Double]
    [vspd,vvxd] = map freturn [vsp,vvx]
    l = V.length vdtd - 1
    [mspds,mvxds] = map (\v -> fromColumns [subVector (g-n) l v | n <- [0..g]]) [vspd,vvxd]
    mx = (konst 1 (l,1)) ! (dropColumns 1 mspds) ! (dropColumns 1 mvxds)
    my = takeColumns 1 mspds
  return (mx,my,vdtd)

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

type PriceData = (String,Double,Double,Double,Double,Double,Double)

type ParsedPriceData = Either String PriceData

readCSV :: FilePath -> IO (V.Vector PriceData)
readCSV f = runSafeT $ runEffect $ PV.runToVectorP $
    (hoist lift (decode HasHeader (PS.withFile f ReadMode (PB.fromHandle))))
    >-> rightP
    >-> PV.toVector

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

mergeEqual :: Ord a => [(a,b)] -> [(a,b)] -> [(a,b,b)]
mergeEqual x [] = []
mergeEqual [] y = []
mergeEqual x@((xi,xv):xs) y@((yi,yv):ys)
  | xi == yi = (xi,xv,yv) : mergeEqual xs ys
  | xi < yi = mergeEqual xs y
  | xi > yi = mergeEqual x ys

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

type Return = Vector Double -> Vector Double

absReturn :: Return
absReturn v = V.zipWith (-) (V.tail v) (V.init v)

pctReturn :: Return
pctReturn v = (V.zipWith (-) (V.tail v) (V.init v))/(V.init v)

cmlAbsReturn :: Return
cmlAbsReturn = V.scanl (+) 0

cmlPctReturn :: Return
cmlPctReturn = V.scanl (\x y -> x * (1+y)) 1

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

readDay :: String -> Day
readDay = read

dayToYear :: Day -> Ratio Integer
dayToYear day =
  let
    (y,_,_) = toGregorian day
    ystart = fromGregorian y 1 1
    yend = fromGregorian y 12 31
    dd = diffDays day ystart
    dy = diffDays yend ystart
  in
    y % 1 + dd % dy

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

rightP :: (Monad m) => Pipe (Either a b) b m r
rightP = for cat $ either (const $ return ()) yield

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

instance MonadPrim m => MonadPrim (SafeT m) where
  type BasePrimMonad (SafeT m) = BasePrimMonad m
  liftPrim = lift . liftPrim

max_stars_repo_name: cmahon/backtests

max_issues_repo_name: cmahon/backtests

max_forks_repo_name: cmahon/backtests


text:
import Data.Complex

import Lib

func c z = z^2 + c

diverge maxIters f x0 = (any (\x -> magnitude x > 4) . take maxIters) (iterate f x0)

points = [a :+ b | a <- [-2.5,-2.495..1], b <- [-1,-0.995..1]]

mset = filter (\x -> not (diverge 1000 (func x) 0)) points

main = plotPoints' "mandelbrot.svg" "Mandelbrot Set" 1.4 (map c2t mset)

max_stars_repo_name: ahaym/dynamics

max_issues_repo_name: ahaym/dynamics

max_forks_repo_name: ahaym/dynamics


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed.Herm
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Hermitian views of packed matrices.
--
module Numeric.LinearAlgebra.Packed.Herm (
    -- * Immutable interface

    -- ** Vector multiplication
    hermMulVector,
    hermMulVectorWithScale,
    addHermMulVectorWithScales,

    -- ** Updates
    hermRank1Update,
    hermRank2Update,

    -- * Mutable interface
    hermCreate,

    -- ** Vector multiplication
    hermMulVectorTo,
    hermMulVectorWithScaleTo,
    addHermMulVectorWithScalesM_,

    -- ** Updates
    hermRank1UpdateM_,
    hermRank2UpdateM_,

    ) where

import Numeric.LinearAlgebra.Packed.Base

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
import Statistics.Distribution
import Statistics.Distribution.Beta

bt :: BetaDistribution
bt = betaDistr a b
       where
         a = 4.0
         b = 2.0

main = do
  putStrLn $ "at least 90 percent: " ++ show (complCumulative bt 0.9)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
-- Copyright (c) 2015-2020 Rudy Matela.
-- Distributed under the 3-Clause BSD licence (see the file LICENSE).
import Test

import System.Exit (exitFailure)
import Data.List (elemIndices)

-- import Test.LeanCheck -- already exported by Test
import Test.LeanCheck.Utils

import Data.Ratio
import Data.Complex
import Data.Int
import Data.Word

main :: IO ()
main  =  do
  max <- getMaxTestsFromArgs 200
  case elemIndices False (tests max) of
    [] -> putStrLn "Tests passed!"
    is -> do putStrLn ("Failed tests:" ++ show is)
             exitFailure

tests :: Int -> [Bool]
tests n =
  [ True

  -- interleave
  , [1,2,3] +| [0,0,0] == [1,0,2,0,3,0]
  , take 3 ([1,2] +| (0:undefined)) == [1,0,2]
  , [0,2..] +| [1,3..] =| n |= [0,1..]

  -- etc
  , tNatPairOrd n
  , tNatTripleOrd n
  , tNatQuadrupleOrd n
  , tNatQuintupleOrd n
  , tNatListOrd n
  , tListsOfNatOrd n
  , listsOf (tiers::[[Nat]]) =| 10 |= tiers

  -- tests!
  , counterExample n (\x y -> x + y /= (x::Int)) == Just ["0", "0"]
  , counterExample n (\x y -> x + y == (x::Int)) == Just ["0", "1"]
  , counterExample n (maybe True (==(0::Int))) == Just ["(Just 1)"]
  , holds n (\x -> x == (x::Int))

  -- For when NaN is in the enumeration (by default, it is not):
  --, fails 100 (\x -> x == (x::Float))  -- NaN != NaN  :-)
  --, counterExample 100 (\x -> x == (x::Float)) == Just ["NaN"]
  , counterExample n (\x y -> x + y == (x::Float))  == Just ["0.0","1.0"]
  , counterExample n (\x y -> x + y == (x::Double)) == Just ["0.0","1.0"]
  , holds          n (\x -> x + 1 /= (x::Int))
  , counterExample n (\x -> x + 1 /= (x::Float))  == Just ["Infinity"]
    || counterExample n (\x -> x + 1 /= (x::Float)) == Just ["inf"] -- bug on Hugs 2006-09?
  , counterExample n (\x -> x + 1 /= (x::Double)) == Just ["Infinity"]
    || counterExample n (\x -> x + 1 /= (x::Float)) == Just ["inf"] -- bug on Hugs 2006-09?
  , allUnique (take n list :: [Float])
  , allUnique (take n list :: [Double])

  , allUnique (take n list :: [Rational])
  , allUnique (take n list :: [Ratio Nat])
  , orderedOn (\r -> numerator r + denominator r) (take n (list :: [Ratio Nat]))
  , orderedOn (\r -> abs (numerator r) + abs(denominator r)) (take n (list :: [Rational]))

  , list == [LT, EQ, GT]
  , orderedOn length (take n (list :: [[Ordering]]))
  , orderedOn length (take n (list :: [[Bool]]))

  , strictlyOrderedOn (\xs -> (sum $ map (+1) xs, xs)) (take n (list :: [[Word]]))

  , tPairEqParams n
  , tTripleEqParams n

  , tProductsIsFilterByLength (tiers :: [[ Nat ]])   10 `all` [1..10]
  , tProductsIsFilterByLength (tiers :: [[ Bool ]])   6 `all` [1..10]
  , tProductsIsFilterByLength (tiers :: [[ [Nat] ]])  6 `all` [1..10]

  , holds n $  (\/)  ==== zipWith' (++) [] [] -:> [[uint2]]
  , holds n $  (\/)  ==== zipWith' (++) [] [] -:> [[bool]]
  , holds n $ (\\//) ==== zipWith' (+|) [] [] -:> [[uint2]]
  , holds n $ (\\//) ==== zipWith' (+|) [] [] -:> [[bool]]

  , holds n $ \x -> x == (x :: Word)
  , holds n $ \x -> x == (x :: Word8)
  , holds n $ \x -> x == (x :: Word16)
  , holds n $ \x -> x == (x :: Word32)
  , holds n $ \x -> x == (x :: Word64)

  , holds n $ \x -> x == (x :: Int)
  , holds n $ \x -> x == (x :: Int8)
  , holds n $ \x -> x == (x :: Int16)
  , holds n $ \x -> x == (x :: Int32)
  , holds n $ \x -> x == (x :: Int64)

  , holds n $ \x -> x == (x :: Complex Double)
  ]

allUnique :: Ord a => [a] -> Bool
allUnique [] = True
allUnique (x:xs) = x `notElem` xs
                && allUnique (filter (< x) xs)
                && allUnique (filter (> x) xs)

-- | 'zipwith\'' works similarly to 'zipWith', but takes neutral elements to
--   operate when one of the lists is exhausted, so, you don't loose elements.
--
-- > zipWith' f z e [x,y] [a,b,c,d] == [f x a, f y b, f z c, f z d]
--
-- > zipWith' f z e [x,y,z] [a] == [f x a, f y e, f z e]
--
-- > zipWith' (+) 0 0 [1,2,3] [1,2,3,4,5,6] == [2,4,6,4,5,6]
zipWith' :: (a->b->c) -> a -> b  -> [a] -> [b] -> [c]
zipWith' _ _  _  []     [] = []
zipWith' f _  zy xs     [] = map (`f` zy) xs
zipWith' f zx _  []     ys = map (f zx) ys
zipWith' f zx zy (x:xs) (y:ys) = f x y : zipWith' f zx zy xs ys

max_stars_repo_name: rudymatela/llcheck

max_issues_repo_name: rudymatela/llcheck

max_forks_repo_name: rudymatela/llcheck


text:
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Statistics
-- Copyright   :  (c) A. V. H. McPhail 2010, 2012, 2014
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Useful statistical functions
--
-----------------------------------------------------------------------------

module Numeric.Statistics (
                           Sample,Samples
                          , covarianceMatrix, correlationCoefficientMatrix
                          , meanList, meanArray, meanMatrix
                          , varianceList, varianceArray, varianceMatrix
                          --
                          , centre, cloglog, corcoeff, cut
                          , ranks, kendall, logit
                          , mahalanobis
                          , mode, moment
                          , ols, percentile, range
                          , run_count
                          , spearman, studentize
                          ) where

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

import Numeric.LinearAlgebra hiding(range)
--import Numeric.LinearAlgebra.Data hiding(range)
--import Numeric.LinearAlgebra.Devel
import Prelude hiding ((<>))
import qualified Data.Array.IArray as I
import qualified Data.List as DL
import qualified Data.Vector.Generic as GV
--import qualified Data.Vector.Storable as SV
import Foreign.Storable
import Numeric.GSL.Statistics
import Numeric.GSL.Sort(sort)

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

type Sample a = Vector a
type Samples a = I.Array Int (Vector a)

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

-- | the covariance matrix
covarianceMatrix :: Samples Double              -- ^ the dimensions of data (each vector being one dimension)
                 -> Matrix Double               -- ^ the symmetric covariance matrix
covarianceMatrix d = let (s,f) = I.bounds d
                      in fromArray2D $ I.array ((s,s),(f,f)) $ concat $ map (\(x,y) -> let c = covariance (d I.! x) (d I.! y) in if x == y then [((x,y),c)] else [((x,y),c),((y,x),c)]) $ filter (\(x,y) -> x <= y) $ I.range ((s,s),(f,f))

-- | the correlation coefficient: (cov x y) / (std x) (std y)
correlationCoefficientMatrix :: Samples Double -> Matrix Double
correlationCoefficientMatrix d = let (s,f) = I.bounds d
                           in fromArray2D $ I.array ((s,s),(f,f)) $ concat $ map (\(x,y) -> let { x' = d I.! x ; y' = d I.! y ; c = (covariance x' y') / ((stddev x') * (stddev y')) } in if x == y then [((x,y),c)] else [((x,y),c),((y,x),c)]) $ filter (\(x,y) -> x <= y) $ I.range ((s,s),(f,f))

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

-- | the mean of a list of vectors
meanList :: (Container Vector a, Num (Vector a), Fractional a) => [Sample a] -> Sample a
meanList []     = error "meanVectors: empty list"
meanList [s]    = s
meanList (s:ss) = let ln = fromIntegral $ length ss + 1
                  in scale (recip ln) $ foldl (+) s ss

-- | the mean of an array of vectors
meanArray :: (Container Vector a, Num (Vector a), Fractional a) => Samples a -> Sample a
meanArray a = meanList $ I.elems a

-- | the mean of a matrix with data series in rows
meanMatrix :: (Container Vector a, Num (Vector a), Fractional a) => Matrix a -> Sample a
meanMatrix a = meanList $ toRows a

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

-- | the variance of a list of vectors
varianceList :: (Container Vector a, Floating (Vector a), Num a, Fractional a) => [Sample a] -> Sample a
varianceList []  = error "varianceList: empty list"
varianceList [s] = konst 0 (size s)
varianceList l   = let mxs = meanList (map (** 2) l)
                       msx = (meanList l) ** 2
                   in mxs - msx

-- | the variance of an array of vectors
varianceArray :: (Container Vector a, Floating (Vector a), Fractional a) => Samples a -> Sample a
varianceArray a = varianceList $ I.elems a

-- | the variance of a matrix with data series in rows
varianceMatrix :: (Container Vector a, Floating (Vector a), Fractional a) => Matrix a -> Sample a
varianceMatrix a = varianceList $ toRows a

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

-- | centre the data to 0: (x - (mean x))
centre :: Vector Double -> Vector Double
centre v = v - (realToFrac (mean v))

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

-- | complementary log-log function
--cloglog :: Vector Double -> Vector Double
cloglog :: Floating a => a -> a
cloglog v = - log (- (log v))

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

-- | corcoeff = covariance x / (std dev x * std dev y)
corcoeff :: Vector Double -> Vector Double -> Double
corcoeff x y = (covariance x y)/((stddev x)*(stddev y))

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

-- | cut numerical data into intervals, data must fall inside the bounds
cut :: Vector Double
    -> Vector Double -- ^ intervals
    -> Vector Int    -- ^ data indexed by bin
cut v c  = let c' = sort c
           in GV.map (\x -> cut_helper 0 x c') v
    where
      cut_helper j x d
          | j >= size d                       = error "Numeric.Statistics: cut: data point not within interval"
          | x >= (d `atIndex` j) && x <= (d `atIndex` (j+1)) = j
          | otherwise                       = cut_helper (j + 1) x d

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

-- | return the rank of each element of the vector
--     multiple identical entries result in the average rank of those entries
--ranks :: Vector Double -> Vector Double
ranks :: (Fractional b, Storable b) => Vector Double -> Vector b
ranks v = let v' = sort v
          in GV.map (\x -> 1 + rank_helper x v') v
              where rank_helper x v' = let is = GV.elemIndices x v'
                                       in (realToFrac (GV.foldl (+) 0 is)) / (fromIntegral $ GV.length is)

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

-- | kendall's rank correlation τ
kendall :: Vector Double -> Vector Double -> Matrix Double
kendall x y = let ln = size x
                  rx = ranks x
                  ry = ranks y
                  r = fromColumns [rx,ry]
                  m = signum $ (kronecker r (asColumn $ konst 1.0 ln)) - (kronecker (asRow $ konst 1.0 ln) r)
                  c = rows m - 1
              in correlationCoefficientMatrix $ I.listArray (0,c) (toColumns m)

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

-- | (logit p) = log(p/(1-p))
--logit :: Vector Double -> Vector Double
logit :: (Floating b, Storable b)
        => Vector b -> Vector b
logit v =  GV.map (\x -> - (log ((1 / x) - 1))) v

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

-- | the Mahalanobis D-square distance between samples
--     columns are components and rows are observations (uses pseudoinverse)
mahalanobis :: Samples Double        -- ^ the data set
            -> Maybe (Sample Double) -- ^ (Just sample) to be measured or use mean when Nothing
            -> Double                -- ^ D^2
mahalanobis x u = let (_,xr) = I.bounds x
                      xl     = I.elems x
                      s'     = pinv $ covarianceMatrix x
                      xu     = case u of
                                 Nothing -> fromList $ map mean xl
                                 Just m  -> m
                      xm     = fromRows $ map ((-) xu) $ toRows $ fromColumns xl
                      --um     = asColumn xu
                      --w      = ((tr' xm) <> xm + (tr' um) <> um)/(fromIntegral $ xr - 1)
                      --w'     = inv w
                  in ((xm <> s' <> (tr' xm)) `atIndex` (0,0))

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

-- | a list of element frequencies
mode :: Vector Double -> [(Double,Integer)]
mode v = let w = sort v
         in DL.sortBy (\(_,n) (_,n') -> compare n' n) $ GV.foldr freqs [] w
            where freqs x []          = [(x,1)]
                  freqs x ((f,n):fns)
                      | f == x         = ((f,n+1):fns)
                      | otherwise     = ((x,1):(f,n):fns)

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

-- | the p'th moment of a vector
moment :: Integral a
       => a             -- ^ moment
       -> Bool          -- ^ calculate central moment
       -> Bool          -- ^ calculate absolute moment
       -> Vector Double -- ^ data
       -> Double
moment p c a v
    | p <= 0     = error "Numeric.Statistics.moment: negative moment requested"
--    | p == 1     = mean v
--    | p == 2     = variance v -- gives sample variance
    | otherwise = let u = if c then centre v else v
                      w = if a then abs u else u
                      x = GV.map (** (fromIntegral p)) w
                  in mean x

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

-- | ordinary least squares estimation for the multivariate model
--   Y = X B + e        rows are observations, columns are elements
--   mean e = 0, cov e = kronecker s I
ols :: (Num (Vector t), Field t)
      => Matrix t         -- ^ X
    -> Matrix t           -- ^ Y
    -> (Matrix t, Matrix t, Matrix t) -- ^ (OLS estimator for B, OLS estimator for s, OLS residuals)
ols x y
    | rows x /= rows y = error "Numeric.Statistics: ols: incorrect matrix dimensions"
    | otherwise       = let (xr,xc) = (rows x,cols x)
                            (yr,yc) = (rows y,cols y)
                            z = (tr' x) <> x
                            r = rank z
                            beta = if r == xc
                                      then (inv z) <> (tr' x) <> y
                                      else (pinv x) <> y
                            rr = y - x <> beta
                            sigma = ((tr' rr) <> rr) / (fromIntegral $ xr - r)
                        in (beta,rr,sigma)

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

-- | compute quantiles in percent
percentile :: Double        -- ^ percentile (0 - 100)
           -> Vector Double -- ^ data
           -> Double        -- ^ result
percentile p d = quantile (0.01*p) d

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

-- | the difference between the maximum and minimum of the input
range :: (Container c e, Num e) => c e -> e
range v = maxElement v - minElement v

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

-- | count the number of runs greater than or equal to @n@ in the data
run_count :: (Num a, Num t, Ord b, Ord a, Container Vector b)
            => a             -- ^ longest run to count
          -> Vector b        -- ^ data
          -> [(a, t)]        -- ^ [(run length,count)]
run_count n v = let w = subVector 1 (size v - 1) v
                    x = GV.foldr run_count' [(1,v `atIndex` 0)] w
                    y = map fst x
                    z = takeWhile (<= n) $  DL.sort y
                in foldr count [] z
    where run_count' m ((c,g):cs)
              | m < g             = ((c+1,m):cs)
              | otherwise         = ((1,m):(c,g):cs)
          run_count' _ [] = error "Lib.Numeric.Statistics.run_count.run_count'"
          count x []           = [(x,1)]
          count x ((yv,yc):ys)
              | x == yv         = ((yv,yc+1):ys)
              | otherwise      = ((x,1):(yv,yc):ys)

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

-- | Spearman's rank correlation coefficient
spearman :: Vector Double -> Vector Double -> Double
spearman x y = corcoeff (ranks x) (ranks y)

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

-- | centre and normalise a vector
studentize :: Vector Double -> Vector Double
studentize x = (centre x)/(fromList $ [stddev x])

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

--table

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

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

max_stars_repo_name: amcphail/hstatistics

max_issues_repo_name: amcphail/hstatistics

max_forks_repo_name: amcphail/hstatistics


text:
{-# LANGUAGE TypeFamilies #-}
module Augmentation.Braid
    (AugBraid (AugBraid)
    ,relations
    ,sChar
    ,s
    ,relations'
    ) where

import Algebra
import Augmentation.DGA
import Braid.Class
import Libs.List

import Data.List
import Data.Either
import Data.Maybe
import Control.Monad
import Data.Functor.Identity

import Numeric.LinearAlgebra

import Debug.Trace

sChar :: Int -> Char
sChar 0 = error "s 0"
sChar k = toEnum $ (2*k)+459

invsChar :: Char -> Int
invsChar c = if fromEnum c > 459 then (fromEnum c - 459) `div` 2 else error "invsChar: c out of bounds"

s :: Int -> Algebra
s = G . E . sChar
{-
sVec :: Int -> Vector R
sVec i = i |> [if j == (i-1) then 1 else 0 | j <- [1..]]
-}
bph :: StdBraid -> Int -> [Int]
bph b@(StdBraid k w) i = nub $ sort $ take k $ iterate (bphh w) i

bphh :: [Int] -> Int -> Int
bphh [] i = i
bphh (w:ws) i
    | w == i = bphh ws (i+1)
    | w == (i-1) = bphh ws (i-1)
    | otherwise = bphh ws i

extractChar :: Algebra -> Maybe Char
extractChar (G (E c)) = Just c
extractChar e = Nothing

sVec :: Int -> Int -> Vector Z
sVec dim i = dim |> [toEnum $ if j == i then 1 else 0 | j <- [1..]]

eqVec :: StdBraid -> Int -> Maybe (Vector Z)
eqVec (StdBraid w ws) i = do
                            { let dim = length ws
                            ; let v = sum $ zipWith (\w x -> if w == i then sVec dim x else if w == (i-1) then negate $ sVec dim x else 0) ws [1..]
                            ; if i > w then Nothing else Just v
                            }

relations' :: StdBraid -> Matrix Z
relations' b@(StdBraid w ws) = maybe (ident $ length ws) id $ do
                                { let basePoints = map head $ nub $ map sort $ map (bph b) [1..w]
                                ; let dim = length ws
                                ; vs <- mapM (eqVec b) $ filter (\i -> not $ i `elem` basePoints) [1..w]
                                ; mat <- if vs == [] then Nothing else Just $ fromColumns vs
                                ; let (q,_) = thinQR $ (fromZ mat :: Matrix R)
                                ; let qqt = q Numeric.LinearAlgebra.<> (tr q)
                                ; n <- (\(a,b) -> if a == b then Just a else Nothing) $ size qqt
                                ; return $ fromColumns $ map (toZ . roundVector) $ toColumns $ 2 * ((ident n) - qqt)
                                }

{-
addVec :: Vector R -> Vector R -> Vector R
addVec v1 v2
    | diff <  0 = addVec v2 v1
    | diff == 0 = v1 + v2
    | otherwise = v1 + (fromList $ toList v2 ++ zers)
    where diff = (size v1) - (size v2)
          zers = take diff repeat 0

eqVec :: StdBraid -> Int -> Maybe (Vector R)
eqVec (StdBraid w ws) i = do
                            { let v = sum $ zipWith (\w x -> if w == i then sVec x else if w == (i-1) then negate $ sVec x else fromInteger 0) ws [1..]
                            ; if i > w then Nothing else Just v
                            }

relations' :: StdBraid -> Maybe (Matrix R)
relations' b@(StdBraid w ws) = do
                                { let vecs = map (eqVec b) [1..w]
                                ; dim <- if [V.length $ head vecs] == (nub $ map V.length vecs) then Just $ V.length $ head vecs else Nothing
                                ; let mat = fmap fromIntegral $ M.transpose $ foldr (\x m -> (colVector x) <|> m) dim vecs
                                ; let zer = zeros (nrows mat) (ncols mat)
                                ; (_,l,p,_) <- luDecomp mat
                                ;

-}
eqh :: StdBraid -> Int -> Maybe (DGA_Map,DGA_Map,DGA_Map,DGA_Map)
eqh (StdBraid _ ws) i = do
                            { let word = zipWith (\w x ->
                                        if w == i
                                            then s x
                                            else if w == (i-1) then recip $ s x
                                                               else 1) ws [1..]
                            ; chars <- mapM (\(x,t) ->
                                    (extractChar x) >>= (\x' -> return (x',t))) $ catMaybes $ map (\(x,t) ->
                                        if x == i then Just (s t,True) else if x == i-1 then Just (s t,False) else Nothing) $ zip ws [1..]
                            ; let (s1,r1) = head chars
                            ; let (s2,r2) = last chars
                            ; let tword = tail word
                            ; let iword = init word
                            ; let prod l = if l == [] then 1 else (head l) * (prod $ tail l)
                            ; let g1 = if r1 then recip else id
                            ; let g2 = if r2 then recip else id
                            ; let m11 = DGA_Map [(s1,g1 $ prod tword)]
                            ; let m12 = DGA_Map [(s1,g1 $ prod $ reverse $ tword)]
                            ; let m21 = DGA_Map [(s2,g2 $ prod iword)]
                            ; let m22 = DGA_Map [(s2,g2 $ prod $ reverse $ iword)]
                            ; return (m11,m12,m21,m22)
                            }

nst :: Int -> (a,a,a,a) -> Maybe a
nst 0 (a,_,_,_) = return a
nst 1 (_,a,_,_) = return a
nst 2 (_,_,a,_) = return a
nst 3 (_,_,_,a) = return a
nst _ _ = Nothing

relations :: StdBraid -> Maybe [DGA_Map]
relations b@(StdBraid w _) = do
                                { let basePoints = map head $ nub $ map sort $ map (bph b) [1..w]
                                ; let solutions = catMaybes $ map (eqh b) $ filter (\x -> not $ x `elem` basePoints) [1..w]
                                --; let toPerm l = foldr (\x xs -> (x `mod` 4) * (4^(length xs)) + xs) l
                                ; let fromPerm m n = if m == 0 then [] else (fromPerm (m-1) ((n - (n `mod` 4)) `div` 4)) ++ [n `mod` 4]
                                ; let perms = map (fromPerm (length solutions)) [0..(4^(length solutions))-1]
                                ; let maps = nub $ catMaybes $ map (\l -> do
                                                    { ms <- zipWithM nst l solutions
                                                    ; let none = DGA_Map []
                                                    ; let left = foldl compose_maps none ms
                                                    ; let right = foldr compose_maps none ms
                                                    ; if left == right && (and $ map (\x -> x == compose_maps x x) ms) then return right else Nothing
                                                    }) perms
                                ; return maps
                                }
{-
eqh :: StdBraid -> Int -> Maybe ((Char,(Algebra,Algebra)),(Char,(Algebra,Algebra)))
eqh b@(StdBraid _ ws) i = do
                            { let word = zipWith (\w x -> if w == i
                                    then s x
                                    else if w == (i-1) then recip $ s x
                                                       else 1) ws [1..]
                            ; c0 <- mapM (\(x,t) -> if x == i then Just (s t,False) else if x == i-1 then Just (s t,True) else Nothing) $ filter (\(x,_) -> (x == i) || (x == i-1)) $ zip ws [1..]
                            ; chars <- mapM (\(x,b) -> (extractChar x) >>= (\x' -> return $ (x',b))) c0
                            ; let (s1,r1) = head chars
                            ; let (s2,r2) = last chars
                            ; let tword = tail word
                            ; let iword = init word
                            ; let pre = ((s1,(product $ tword, product $ reverse $ tword)),(s2,(product $ iword,product $ reverse $ iword)))
                            ; let g1 = if r1 then recip else id
                            ; let g2 = if r2 then recip else id
                            ; let f ((c1,(e11,e12)),(c2,(e21,e22))) = ((c1,(g1 e11,g1 e12)),(c2,(g2 e21, g2 e22)))
                            ; return $ f pre
                            }

relations :: StdBraid -> Maybe [DGA_Map]
relations b@(StdBraid w _) = do
                                { let basePoints = map head $ nub $ map sort $ map (bph b) [1..w]
                                ; solutions <- mapM (eqh b) $ filter (\x -> not $ x `elem` basePoints) [1..w]
                                ; let solution = catMaybes $ map (\((c1,(a11,a12)),(c2,(a21,a22))) -> if c2 == c1 then Nothing else Just [(c1,a11),(c1,a12),(c2,a21),(c2,a22)]) solutions
                                ; let splited = map (\x -> compose_maps x x) $ map (\x -> DGA_Map $ zipWith (\l x' -> l !! x') solution x) $ permutations [0..3]
                                ; let fixed = nub $ filter (\x -> x == (compose_maps x x)) splited
                                ; return fixed
                                }
-}
buildMaph :: Maybe [(Char,Algebra)] -> Maybe [(Char,Algebra)] -> Maybe [(Char,Algebra)]
buildMaph m m'
    | m == m' = m
    | otherwise = buildMaph (m >>= (\l -> Just $ map (\(c,x) -> (c,applyDGAMap (DGA_Map l) x)) l)) m

footprinth :: Int -> [Either (Int,Char) Int] -> [(FreeGroup,Int)]
footprinth _k [] = []
footprinth k ((Right i):xs) = (E $ toEnum k,i):(footprinth (k+1) xs)
footprinth k ((Left (i,c)):xs) = [(E c,i-1),(invert $ E c,i+1)] ++ (footprinth (k+1) xs)

iscrossh :: [Either (Int,Char) Int] -> Int -> Bool
iscrossh [] _ = False
iscrossh ((Left _):_) 0 = False
iscrossh ((Left _):_) 1 = False
iscrossh ((Right _):_) 0 = True
iscrossh ((Right _):cs) x = iscrossh cs (x-1)
iscrossh ((Left _):cs) x = iscrossh cs (x-2)

data AugBraid = AugBraid Int [Either (Int,Char) Int] -- Either element of H1(L), and its position or an integer representing the corresponding braid group element
instance Braid AugBraid where
    type M AugBraid = Either (Int,Char)
    get_word (AugBraid _ w) = w
    algebra_footprint (AugBraid _ w) = map (\(g,i) -> (G g,i)) $ footprinth 65 w
    toStdBraid (AugBraid w ws) = StdBraid w (rights ws)
    fromStdBraid (StdBraid w ws) = AugBraid w (map Right ws)
    isCross (AugBraid _ ws) x = iscrossh ws x
    cross_art b row (Right s) = cross_art (toStdBraid b) row (Identity s)
    cross_art _b row (Left (s,s')) = let s1 = head $ show $ E s'
                                         s2 = show $ invert $ E s'
                                      in if (row - (abs $ s-1)*3) `elem` [0..3]
                        then Just $ (case ((row-(s-1)*3) `mod` 4) of 0 -> "---" ++ [s1] ++ "----"
                                                                     1 -> "   " ++  " " ++ "    "
                                                                     2 -> "   " ++  " " ++ "    "
                                                                     3 -> "--" ++ s2 ++ "---")
                        else Nothing
instance Eq AugBraid where
    (==) = equal
instance Show AugBraid where
    show = showBraid

max_stars_repo_name: Creatorri/Legendrian-Knots-UROP

max_issues_repo_name: Creatorri/Legendrian-Knots-UROP

max_forks_repo_name: Creatorri/Legendrian-Knots-UROP


text:
module LinReg (

  lsr, bulkPDs

  ) where

import Poly
import Data.List
import Data.Monoid
import Control.Arrow
import Data.Packed.Matrix
import Numeric.LinearAlgebra.LAPACK

type Data = [(Double,Double)]

-- Data Acquisition - Very unsafe

getData :: PartialDerivitives -> ([[Double]], [Double])
getData = (map tail &&& map (negate . head)) . fun
  where
    fun                         = map (map (getNums . unSP) . unSS) . unPD
    getNums                     = foldl' getNum 1
    getNum x (SimpleNumber   n) = n * x
    getNum x (SimpleConstant _) = x

-- Concat Partialderivitives over data-points

bulkPDs :: Int -> Data -> PartialDerivitives
bulkPDs ord dat = pard
  where
    poly = polynomial_of_order               ord
    dsqr = delta_squared_term                poly
    ssss = mconcat $ map (point_result dsqr) dat
    pard = partial_derivatives ord           ssss

point_result :: Expression -> (Double,Double) -> SimpleSum
point_result e (x,y) = normal (Env [("x", x), ("y", y)]) e

lsr :: Int -> Data -> [Double]
lsr order points = concat $ toLists $ linearSolveR coefs values
  where
    linear     = bulkPDs (order - 1) points
    (lhs, rhs) = getData linear
    coefs      = (order >< order) (concat lhs)
    values     = (order >< 1    )  rhs

max_stars_repo_name: sordina/Simple-LSR

max_issues_repo_name: sordina/Simple-LSR

max_forks_repo_name: sordina/Simple-LSR


text:
module NeuralNet.Example (
  Example,
  ExampleSet (),
  createExampleSet,
  exampleSetX,
  exampleSetY,
  exampleSetM,
  exampleSetN,
  splitExampleSet
) where

import Numeric.LinearAlgebra
import Data.List (nub)

type Example = ([Double], Double)

data ExampleSet = ExampleSet [Example] (Matrix Double) (Matrix Double)
  deriving (Show, Eq)

createExampleSet :: [Example] -> ExampleSet
createExampleSet e
  | numUniqueExampleSizes > 1  = error ("All examples must have same size (given " ++ show uniqueExampleSizes ++ ")")
  | numUniqueExampleSizes == 0 = error "No examples provided"
  | n == 0                     = error "Empty example provided"
  | otherwise                  = ExampleSet e x y
    where
      exampleSizes = map (length . fst) e
      uniqueExampleSizes = nub exampleSizes
      numUniqueExampleSizes = length uniqueExampleSizes
      n = head exampleSizes
      x = tr (fromLists (map fst e))
      y = fromLists [map snd e]

exampleSetX :: ExampleSet -> Matrix Double
exampleSetX (ExampleSet _ x _) = x

exampleSetY :: ExampleSet -> Matrix Double
exampleSetY (ExampleSet _ _ y) = y

exampleSetN :: ExampleSet -> Int
exampleSetN = rows . exampleSetX

exampleSetM :: ExampleSet -> Int
exampleSetM (ExampleSet e _ _) = length e

splitExampleSet :: Double -> ExampleSet -> (ExampleSet, ExampleSet)
splitExampleSet ratio (ExampleSet e _ _) = (createExampleSet leftExamples, createExampleSet rightExamples)
  where
    (leftExamples, rightExamples) = splitAt (round (ratio * fromIntegral (length e))) e


text:
{-# LANGUAGE CPP                   #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Trustworthy
           , DefaultSignatures     #-}
#define USE_GHC_GENERICS
#endif
{-# LANGUAGE RankNTypes
           , TypeFamilies
           , KindSignatures
           , DeriveGeneric
           , FlexibleInstances
           , FlexibleContexts
           , BangPatterns
           , MultiParamTypeClasses #-}
{-# OPTIONS_GHC -fno-warn-orphans  #-}

-------------------------------------------------------------------------------
-- |
-- Module      : OGLS.Engine.Math.Vectors
-- description : 3D and 4D vectors arithmetic
-- Copyright   : (c) Adam, 2017
-- License     : MIT
-- Maintainer  : awkure@protonmail.ch
-- Stability   : unstable
-- Portability : POSIX
-------------------------------------------------------------------------------

module OGLS.Engine.Math.Vectors
    ( (^$+), (^$-), (^$*)
    , dot, step
    ) where

#ifdef USE_GHC_GENERICS
import GHC.Generics               ( Generic (..) )
#endif
import Control.Applicative        ( liftA2       )
import Control.Monad              ( liftM2       )
import Control.Parallel           ( pseq         )
import Data.Complex               ( Complex (..) )
import Data.Orphans               (              )
import Control.Lens hiding ( (<.>) )

import OGLS.Engine.Math.Instances

infixl 6 ^$+, ^$-
infixl 7 ^$*

default ()

class Vector a v where
    -- not yet implemented

(^$+) :: forall (t :: * -> *) a. (Applicative t, Num a) => t a -> t a -> t a
(^$-) :: forall (t :: * -> *) a. (Applicative t, Num a) => t a -> t a -> t a
(^$*) :: forall (t :: * -> *) a. (Applicative t, Num a) => t a -> t a -> t a

newtype E t = E { el :: forall x. Lens' (t x) x }

data Vec3 a = Vec3 { v3x :: {-# UNPACK #-} !a
                   , v3y :: {-# UNPACK #-} !a
                   , v3z :: {-# UNPACK #-} !a
                   } deriving ( Eq, Ord, Show, Generic )

instance Functor Vec3 where
    fmap f (Vec3 x y z) = Vec3 (f x) (f y) (f z)

instance Applicative Vec3 where
    pure x = Vec3 x x x
    Vec3 f g h <*> Vec3 x y z = Vec3 (f x) (g y) (h z)

data Vec4 a = Vec4 { v4x :: {-# UNPACK #-} !a
                   , v4y :: {-# UNPACK #-} !a
                   , v4z :: {-# UNPACK #-} !a
                   , v4w :: {-# UNPACK #-} !a
                   } deriving ( Eq, Ord, Show )

(^$+) = liftA2 (+)
{-# INLINE (^$+) #-}
(^$-) = liftA2 (-)
{-# INLINE (^$-) #-}
(^$*) = liftA2 (*)
{-# INLINE (^$*) #-}

-- | Left scalar product vector
(*^^) :: forall (f :: * -> *) a. (Functor f, Num a) => a -> f a -> f a
(*^^) a = fmap (a*)
{-# INLINE (*^^) #-}

-- | Right scalar product vector
(^^*) :: forall (f :: * -> *) a. (Functor f, Num a) => f a -> a -> f a
f ^^* a = fmap (*a) f
{-# INLINE (^^*) #-}

{- Needs Additive class and instances
basis :: (Traversable t, Num a) => [t a]
basis = basisFor (zero :: Additive v => v Int)

basisFor :: (Traversable t, Num a) => t b -> [t a]
basisFor = \t ->
   ifoldMapOf traversed ?? t $ \i _ ->
     return $ iover  traversed ?? t $ \j _ ->
         if i == j then 1 else 0
{-# INLINABLE basisFor #-}
-}

-- | Linear interpolation
lerp :: forall (f :: * -> *) a. (Functor f, Applicative f, Num a) => a -> f a -> f a -> f a
lerp alpha u v = alpha *^^ u ^$+ (1 - alpha) *^^ v
{-# INLINE lerp #-}

-- | Dot product
dot :: forall (t :: * -> *) a. (Applicative t, Foldable t, Num a) => t a -> t a -> a
dot v1 v2 = sum $ v1 ^$* v2

step :: forall (t :: * -> *) a. (Applicative t, Num a, Ord a) => t a -> t a -> t a
step = liftA2 (\a b -> if b < a then 0 else 1)

#ifdef HLINT
{-# ANN module "HLint: ignore Redundant lambda" #-}
#endif

max_stars_repo_name: awkure/ogls

max_issues_repo_name: awkure/ogls

max_forks_repo_name: awkure/ogls


text:
module Fractal.Math.Kernel(
  evaluateFractal
)where

import Data.Complex

type Point = Complex Double
type Limit = Double
type Iteration = Complex Double
type Magnitude = Double
type MaxIter = Int
type Iter = Int

evaluateFractal :: (Point -> Iteration -> Iteration) -> Limit -> MaxIter -> Point -> Iter
evaluateFractal f l m c = res c (0.0 :+ 0.0) 0
  where res c z i             = if i == m then i else rest (f c z) i
        rest val i             = if checkAcceptance l (magnitude val) then i else res c val (i + 1)

checkAcceptance :: Limit -> Double -> Bool
checkAcceptance l val
  | val >= l    = True
  | val == 0.0  = True
  | otherwise   = False

max_stars_repo_name: sigmaticsMUC/FractalHaskell

max_issues_repo_name: sigmaticsMUC/FractalHaskell

max_forks_repo_name: sigmaticsMUC/FractalHaskell


text:
{-# LANGUAGE FlexibleContexts       #-}
{-# LANGUAGE FlexibleInstances      #-}
{-# LANGUAGE MultiParamTypeClasses  #-}
{-# LANGUAGE QuantifiedConstraints  #-}
{-# LANGUAGE RecordWildCards        #-}
{-# LANGUAGE ScopedTypeVariables    #-}
{-# LANGUAGE UndecidableInstances   #-}
module Q.Stochastic.Discretize
        where

import           Data.Functor
import           Numeric.LinearAlgebra
import           Q.Stochastic.Process
-- |Euler discretization of stochastic processes
newtype Euler = Euler { eDt :: Double }
        deriving stock (Show, Eq)

-- | Euler end-point discretization of stochastic processes
newtype EndEuler = EndEuler { eeDt :: Double }
        deriving stock (Show, Eq)

instance Discretize Euler Double where
  dDrift p Euler{..} s0 = pDrift p s0 <&> (* eDt)
  dDiff  p Euler{..} b  = (pDiff p b) <&> (* (sqrt eDt))
  dDt    _ Euler{..} _  = eDt

instance Discretize Euler (Vector Double) where
  dDrift p Euler{..} s0 = pDrift p s0 <&> (scale eDt)
  dDiff  p Euler{..} b = (pDiff p b) <&> (scale (sqrt eDt))
  dDt    _ Euler{..} _  = eDt

instance Discretize EndEuler Double where
  dDrift p EndEuler{..} s0@(t0, x0) = pDrift p (t0 + eeDt, x0) <&> (* eeDt)
  dDiff  p EndEuler{..}  s0@(t0, x0) =  pDiff  p (t0 + eeDt, x0) <&> (* (sqrt eeDt))
  dDt    _ e _   = eeDt e

max_stars_repo_name: ghais/HQu

max_issues_repo_name: ghais/HQu

max_forks_repo_name: ghais/HQu


text:
{- Assignment 1
 - Name: Omar Alkersh
 - Date: 22/09/2018
 -}
module Assign_1 where
import           Data.Complex

macid :: String
macid = "alkersho"
{-
 - Throughout this file I will be referring to https://proofwiki.org/wiki/Cardano%27s_Formula
 - and the variables listed within it.
 - Equation in the form of: ax^3 + bx^2 + cx + d = 0
-}
{- -----------------------------------------------------------------
 - cubicQ
 - -----------------------------------------------------------------
 - Description: Finds Q, refer to website
 -}
cubicQ :: Double -> Double -> Double -> Double
cubicQ a b c = ((3*a*c)-b^2)/(9*a^2)

{- -----------------------------------------------------------------
 - cubicR
 - -----------------------------------------------------------------
 - Description: Finds R using
 -}
cubicR :: Double -> Double -> Double -> Double -> Double
cubicR a b c d = ((9*a*b*c)-(27*a^2*d)-(2*b^3))/(54*a^3)

{- -----------------------------------------------------------------
 - cubicDisc
 - -----------------------------------------------------------------
 - Description: Finds the discriminant
 - when the disc > 0, there is one real solution
 - when the disc = 0, there are two real solutions
 - when the disc < 0, there are three real solutions
 -}
cubicDisc :: Double -> Double -> Double
cubicDisc q r = q^3 + r^2

--Returns the cube root of a number
cubeRoot :: Double -> Double
cubeRoot x = x ** (1/3)

{- -----------------------------------------------------------------
 - cubicS
 - -----------------------------------------------------------------
 - Description: Finds S using Q and R
 -}
cubicS :: Double -> Double -> Double
cubicS q r = cubeRoot(r+sqrt(q^3 + r^2))

{- -----------------------------------------------------------------
 - cubicT
 - -----------------------------------------------------------------
 - Description: Finds T using Q and R
 -}
cubicT :: Double -> Double -> Double
cubicT q r = cubeRoot(r-sqrt(q^3 + r^2))

{- -----------------------------------------------------------------
 - cubicRealSolutions
 - -----------------------------------------------------------------
 - Description: Finds the roots/solutions of the cubic equaition
 - when the disc > 0, there is one real solution
 - when the disc = 0, there are two real solutions
 - when the disc < 0, there are three real solutions
 -}
cubicRealSolutions :: Double -> Double -> Double -> Double -> [Double]
cubicRealSolutions a b c d =
  let
    q = cubicQ a b c
    r = cubicR a b c d
    s = cubicS q r
    t = cubicT q r
    x1 = s + t - (b / (3*a))
    x2 = -((s+t)/2)-(b/3*a)+((sqrt 3 / 2)*(s-t))
    tol = 1e-5
    in if cubicDisc q r > 0
     then [x1]
     else if abs (cubicDisc q r) == 0.0
       then [x1, x2, x2]
       else []
{- -----------------------------------------------------------------
 - Test Cases
 - -----------------------------------------------------------------
 -}

--test returns True if test succeds
--Formulae in question:
-- x^3 + 3x^2 + 3x + 1 = (x+1)^3
-- x^3 - 6x^2 + 9x + 0 = x(x-3)^2
test :: (Double, Double) -> (Double, Double) -> (Double, Double) -> Bool
test xs1 xs2 xs3 = let
 a = fst xs1 * fst xs2 * fst xs3
 b = fst xs1 * fst xs2 * snd xs3 + fst xs1 * snd xs2 * fst xs3 + snd xs1 * fst xs2 * fst xs3
 c = fst xs1 * snd xs2 * snd xs3 + snd xs1 * fst xs2 * snd xs3 + snd xs1 * snd xs2 * fst xs3
 d = snd xs1 * snd xs2 * snd xs3
 x1 = -(snd xs1 / fst xs1)
 x2 = -(snd xs2 / fst xs2)
 x3 = -(snd xs3 / fst xs3)
 in cubicRealSolutions a b c d == [x1, x2, x2] || cubicRealSolutions a b c d == [] || cubicRealSolutions a b c d == [x1]