text:
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ViewPatterns          #-}

-- | TODO
module SatOpt.Render ( Conf (..)
                     , runDisp
                     , KMState --(..)
                     , RGBTrips
                     ) where

import           Control.Concurrent.MVar
import           Control.Monad                 (when)
import           Data.Maybe                    (isJust)
import           Data.Packed.Vector            (buildVector)
import           Foreign                       (withArray)
import           Graphics.UI.GLUT
import           Numeric.LinearAlgebra.HMatrix hiding (reshape)
import           System.Exit                   (exitSuccess)

type Mℝ = Numeric.LinearAlgebra.HMatrix.Matrix ℝ

type RGBTrips = (Mℝ, Mℝ, Mℝ)

type KMState = String

data Conf a = Conf { keyBinds :: KMState -> (Conf a) -> IO (Conf a)
                   , state    :: a
                   , evolve   :: a -> a
                   , render   :: a -> RGBTrips
                   , canvas   :: (Int, Int)
                   }

imageSize :: Conf a -> TextureSize2D
imageSize (canvas -> (x, y)) = TextureSize2D (fromIntegral x) (fromIntegral y)

(//) :: (Integral a, Fractional b) => a -> a -> b
a // b = fromIntegral a / fromIntegral b

toColor :: Double -> Double -> Double -> Color4 GLubyte
toColor r g b = Color4 (tc r) (tc g) (tc b) 255
  where
    tc = round . (*250)

mapVecs3 :: (Container Vector a, Indexable (Vector a) a, Element b)
            => (a -> a -> a -> b)
            -> Vector a
            -> Vector a
            -> Vector a
            -> Vector b
mapVecs3 f x y z = buildVector (size x) (\i -> f (x ! i) (y ! i) (z ! i))

instance Element GLubyte
instance Element (Color4 GLubyte)

flattenImage :: RGBTrips -> Vector (Color4 GLubyte)
flattenImage (flatten -> r, flatten -> g, flatten -> b)
  = mapVecs3 toColor r g b

withImage :: Conf a -> (PixelData (Color4 GLubyte) -> IO ()) -> IO ()
withImage c act = withArray arr $ act . PixelData RGBA UnsignedByte
  where
    arr = toList $ flattenImage $ render c (state c)

myInit :: Conf a -> IO (Maybe TextureObject)
myInit c = do
   clearColor $= Color4 0 0 0 0
   shadeModel $= Flat
   depthFunc $= Just Less
   rowAlignment Unpack $= 1

   exts <- get glExtensions
   mbTexName <- if "GL_EXT_texture_object" `elem` exts
                   then fmap Just genObjectName
                   else return Nothing
   when (isJust mbTexName) $ textureBinding Texture2D $= mbTexName

   textureWrapMode Texture2D S $= (Repeated, Repeat)
   textureWrapMode Texture2D T $= (Repeated, Repeat)
   textureFilter Texture2D $= ((Nearest, Nothing), Nearest)
   withImage c $ texImage2D Texture2D NoProxy 0 RGBA' (imageSize c) 0
   texture Texture2D $= Enabled
   return mbTexName

evolveC :: Conf a -> Conf a
evolveC c = c { state = evolve c $ state c }

display :: MVar (Conf a) -> Maybe TextureObject -> DisplayCallback
display mv mbTexName = do
  cstate <- readMVar mv
  clear [ ColorBuffer, DepthBuffer ]
  textureFunction $= Decal
  when (isJust mbTexName) $ textureBinding Texture2D $= mbTexName
  let texCoord2f = texCoord :: TexCoord2 GLfloat -> IO ()
      vertex3f   = vertex   :: Vertex3   GLfloat -> IO ()
  renderPrimitive Quads $ do
    texCoord2f (TexCoord2 0 0); vertex3f (Vertex3 (-1.0) (-1.0) 0.0 )
    texCoord2f (TexCoord2 0 1); vertex3f (Vertex3 1.0 (-1.0) 0.0 )
    texCoord2f (TexCoord2 1 1); vertex3f (Vertex3 1.0 1.0 0.0 )
    texCoord2f (TexCoord2 1 0); vertex3f (Vertex3 (-1.0) 1.0 0.0 )
  withImage cstate $ texImage2D Texture2D NoProxy 0 RGBA' (imageSize cstate) 0
  swapBuffers

reshape :: ReshapeCallback
reshape size@(Size w h) = do
  viewport $= (Position 0 0, size)
  matrixMode $= Projection
  loadIdentity
  perspective 60 (fromIntegral w / fromIntegral h) 1 1
  matrixMode $= Modelview 0
  loadIdentity
  translate (Vector3 0 0 (-1 :: GLfloat))

kmEncode :: Key -> KeyState -> Modifiers -> Position -> KMState
kmEncode (Char k)        Down _ _ = [k]
kmEncode (SpecialKey k)  Down _ _ = show k
kmEncode (MouseButton m) Down _ _ = show m
kmEncode _ _ _ _ = ""

keyboard :: MVar (Conf a) -> KeyboardMouseCallback
keyboard _ (Char '\27') Down _ _ = exitSuccess
keyboard mv k ks m p = do
  c <- readMVar mv
  let e = kmEncode k ks m p
  modifyMVar_ mv (keyBinds c e)

idle :: MVar (Conf a) -> IdleCallback
idle mv = do
  modifyMVar_ mv (return . evolveC)
  postRedisplay Nothing

runDisp :: Conf a -> IO ()
runDisp c = do
  cstate <- newMVar c
  (progName, _args) <- getArgsAndInitialize
  initialDisplayMode    $= [ DoubleBuffered, RGBMode, WithDepthBuffer ]
  initialWindowSize     $= Size 250 250
  initialWindowPosition $= Position 100 100
  _ <- createWindow progName
  tex <- myInit c
  displayCallback       $= display cstate tex
  reshapeCallback       $= Just reshape
  keyboardMouseCallback $= Just (keyboard cstate)
  idleCallback          $= Just (idle cstate)
  mainLoop

{-

module Main where

    import Graphics.GPipe
    import qualified Data.Vec as Vec
    import Data.Vec.Nat
    import Data.Vec.LinAlg.Transform3D
    import Data.Monoid
    import Data.IORef
    import Foreign.C.Types
    import Foreign.Marshal.Array
    import Graphics.UI.GLUT( Window,
                             mainLoop,
                             postRedisplay,
                             idleCallback,
                             getArgsAndInitialize,
                             ($=))

    triVerts = [ (-0.8):.0.8:.0.0:.(),
             0.8:.0.8:.0.0:.(),
             0.0:.(-0.8):.0.0:.() ]
    uvCoords = [ 0.0:.0.0:.(),
                 1.0:.0.0:.(),
                 0.5:.1.0:.() ]

    triangle :: PrimitiveStream Triangle (Vec3 (Vertex Float),Vec2 (Vertex Float))
    triangle = toGPUStream TriangleList $ zip triVerts uvCoords

    -- This implements the vertex shader
    procTriangle :: Vec2 Int -> PrimitiveStream Triangle (Vec4 (Vertex Float), Vec2 (Vertex Float))
    procTriangle size = fmap (projTriangle size) triangle

    projTriangle :: Vec2 Int -> (Vec3 (Vertex Float),Vec2 (Vertex Float)) -> (Vec4 (Vertex Float), Vec2 (Vertex Float))
    projTriangle size (pos,uv) = (orthoProj `multmv` homPos, uv)
        where homPos = homPoint pos :: Vec4 (Vertex Float)
              orthoProj = toGPU $ orthogonal (-10) 10  (2:.2:.())

    -- This implements the fragment shader
    rastTriangle :: Texture2D RGBFormat -> Vec2 Int -> FragmentStream (Color RGBFormat (Fragment Float))
    rastTriangle tex size = fmap (applyTex tex) $ rasterizeFrontAndBack $ procTriangle size

    applyTex tex (front,uv) = sample (Sampler Linear Wrap) tex uv

    triangleFrame :: Texture2D RGBFormat -> Vec2 Int -> FrameBuffer RGBFormat () ()
    triangleFrame tex size = draw (rastTriangle tex size) clear
        where
          draw  = paintColor NoBlending (RGB $ Vec.vec True)
          clear = newFrameBufferColor (RGB (0.1:.0.3:.0.6:.()))

    main :: IO ()
    main = do
      getArgsAndInitialize
      putStrLn "loading texture..."
      tex <- loadTexture demonHead
      putStrLn "creating window..."
      newWindow "Green Triangle"
           (100:.100:.())
           (800:.600:.())
           (renderFrame tex)
           initWindow
      putStrLn "entering mainloop..."
      mainLoop

    renderFrame :: Texture2D RGBFormat -> Vec2 Int -> IO (FrameBuffer RGBFormat () ())
    renderFrame tex size = do
              return $ triangleFrame tex size

    initWindow :: Window -> IO ()
    initWindow win = idleCallback $= Nothing --Just (postRedisplay (Just win))

    loadTexture' ptr = newTexture (PerComp3 UnsignedByteFormat) RGB8 (128:.128:.()) [ptr]
    loadTexture = flip withArray loadTexture'

    demonHead :: [CUChar]

-}

max_stars_repo_name: taktoa/SatOpt

max_issues_repo_name: taktoa/SatOpt

max_forks_repo_name: taktoa/SatOpt


text:
{-# LANGUAGE OverloadedStrings
  , FlexibleInstances
  , GeneralizedNewtypeDeriving
  , MultiParamTypeClasses
  , InstanceSigs
  , QuasiQuotes #-}
module Main where

import Data.List

import Control.Monad.IO.Class
import Control.Monad.Trans
import Control.Monad.Reader

import Data.Maybe
import qualified Data.Vector as V
import Data.Time

import Types
import Persistence

import Colog
import Colog.Message
import Colog.Actions

import Graphics.Rendering.Chart.Easy hiding ((<.>))
import Graphics.Rendering.Chart.Backend.Diagrams (toFile)

import System.Environment
import System.FilePath.Posix
import Statistics.Quantile

import Text.InterpolatedString.Perl6 (qc)

data Env m = Env { _persistence :: Persistence
                 , _logAction :: !(LogAction m Message) }

newtype App e = App { unapp :: ReaderT (Env App) IO e }
    deriving (Monad, Functor, Applicative, MonadIO, MonadReader (Env App))

instance HasPersistence (Env m) where
    getPersistence = _persistence

instance HasLog (Env m) Message m where
    getLogAction = _logAction
    setLogAction newLogAction env = env { _logAction = newLogAction }

indexHtml diagram = [qc|

        <title>SMHI Forecast error</title>

|]

withApp :: Env App -> App () -> IO ()
withApp env f = runReaderT (unapp f) env

-- environment variable with default value in io monad
fromEnv :: String -> String -> IO String
fromEnv d = fmap (fromMaybe d) . lookupEnv

type BoxValues = (Double, Double, Double, Double, Double)

lineStyle :: Double -> Colour Double -> LineStyle
lineStyle n colour = line_width .~ n
                   $ line_color .~ opaque colour
                   $ def

box label color vals = liftEC $ do
    plot_candle_line_style  .= lineStyle 1 color
    plot_candle_fill .= False
    plot_candle_rise_fill_style .= solidFillStyle (opaque color)
    plot_candle_fall_fill_style .= solidFillStyle (opaque color)
    plot_candle_tick_length .= 2
    plot_candle_width .= 8
    plot_candle_centre .= 8
    plot_candle_values .= [ Candle ind mn first med third mx | (ind,(mn,first,med,third,mx)) <- vals]
    plot_candle_title .= label

grpBox :: [HourOffset] -> (Double, (Double, Double, Double, Double, Double))
grpBox os = (k, (mn, q1, q2, q3, mx))
    where k = fromIntegral $ hoursAgo (head os)
          mn = quantile s 1 10 temperatures
          mx = quantile s 9 10 temperatures
          q1 = quantile s 1 4 temperatures
          q2 = quantile s 2 4 temperatures
          q3 = quantile s 3 4 temperatures
          temperatures = V.fromList (map temperatureE os)

every n xs = case drop (n-1) xs of
              (y:ys) -> y : every n ys
              [] -> []

errorsByHour = groupBy (\a b -> hoursAgo a == hoursAgo b) . sortOn hoursAgo

main :: IO ()
main = do
    time <- show <$> getCurrentTime
    outdir <- fromEnv "./" "OUTDIR"
    database <- fromEnv "data.db" "DATABASE"
    image <- fromEnv (outdir </> time <.> "svg") "DIAGRAM"
    index <- fromEnv (outdir </> "index.html") "HTML"

    env <- Env <$> initPersistence database
               <*> pure richMessageAction

    withApp env $ do
        samples <- getHourOffsets
        let grps = errorsByHour samples
            positives = map (filter (\a -> temperatureE a > 0)) grps
            negatives = map (filter (\a -> temperatureE a < 0)) grps
            positivesBoxes = map grpBox $ filter (\grp -> length grp > 2) positives
            negativesBoxes = map grpBox $ filter (\grp -> length grp > 2) negatives

        mae <- getMAE
        liftIO $ toFile def image $ do
            layout_title .= "SMHI Forecast Accuracy (Generated " ++ show time ++")"
            layout_y_axis . laxis_title   .= "Temperature error (C)"
            layout_x_axis . laxis_title   .= "Forecast age (h)"
            layout_x_axis . laxis_reverse .= True

            plot (box "high forecast" blue $ every 10 negativesBoxes)
            plot (box "low forecast"  red  $ every 10 positivesBoxes)

        liftIO $ writeFile index (indexHtml image)

ext: hs

max_stars_repo_name: kaffepanna/smhi

max_issues_repo_name: kaffepanna/smhi

max_forks_repo_name: kaffepanna/smhi


text:
import Data.Complex (cis, phase)
import Data.List.Split (splitOn)
import Text.Printf (printf)

timeToRadians :: String -> Float
timeToRadians time =
        let hours:minutes:seconds:_ = splitOn ":" time
            s = fromIntegral (read seconds :: Int)
            m = fromIntegral (read minutes :: Int)
            h = fromIntegral (read hours :: Int)
        in  (2*pi)*(h+ (m + s/60.0 )/60.0 )/24.0

radiansToTime :: Float -> String
radiansToTime  r =
        let tau = pi*2
            (_,fDay) = properFraction (r / tau) :: (Int, Float)
            fDayPositive = if fDay < 0 then 1.0+fDay else fDay
            (hours, fHours) = properFraction $ 24.0 * fDayPositive
            (minutes, fMinutes) = properFraction $ 60.0 * fHours
            seconds = 60.0 * fMinutes
        in printf "%0d" (hours::Int) ++ ":" ++ printf "%0d" (minutes::Int) ++ ":" ++ printf "%0.0f" (seconds::Float)

meanAngle :: [Float] -> Float
meanAngle = phase . sum . map cis

main :: IO ()
main = putStrLn $ radiansToTime $ meanAngle $ map timeToRadians ["23:00:17", "23:40:20", "00:12:45", "00:17:19"]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Qubism.StateVecSpec (genStateVec, spec) where

import GHC.TypeLits
import Data.Singletons

import           Test.Hspec
import           Test.QuickCheck
import           Test.QuickCheck.Monadic
import           Control.Monad.Random
import           Control.Monad.Trans.State
import qualified Numeric.LinearAlgebra as LA
import           Data.Complex

import Qubism.AlgebraTests
import Qubism.StateVec

genC :: Gen (Complex Double)
genC = do
  a <- choose (-1, 1)
  b <- choose (-1, 1)
  pure $ a :+ b

genStateVec :: forall n . KnownNat n => Gen (StateVec n)
genStateVec = normalize . UnsafeMkStateVec . LA.fromList <$> vectorOf l genC
  where l = (2 ^) . fromIntegral . fromSing $ (sing :: Sing n)

-- | I know, it's an orphan instance, but it's not relevant outside of this
-- test suite, and never should be.
instance KnownNat n => Arbitrary (StateVec n) where
  arbitrary = genStateVec

propIdempotent
  :: (Eq a, Show a, Eq b, Show b)
  => StateT a (Rand StdGen) b
  -> a
  -> PropertyM IO Expectation
propIdempotent st a = do
  g <- lift getStdGen
  let one = runStateT st a `evalRand` g
      two = runStateT (st >> st) a `evalRand` g
  pure $ two `shouldBe` one

spec :: Spec
spec = do
  describe "Quantum Registers" $ do
    describe "form a vector space"  $ isVectorSpace  (T :: T (StateVec 1))
    describe "form a hilbert space" $ isHilbertSpace (T :: T (StateVec 1))
    describe "measure" $
      it "is idempotent"
        $ property
        $ monadicIO
        $ forAllM (genStateVec :: Gen (StateVec 1))
        $ propIdempotent measure
    describe "measureQubit" $
      it "is idempotent"
        $ property
        $ monadicIO
        $ forAllM (genStateVec :: Gen (StateVec 1))
        $ propIdempotent (measureQubit 0)

max_stars_repo_name: qubitrot/qubism

max_issues_repo_name: qubitrot/qubism

max_forks_repo_name: qubitrot/qubism


text:
module Numeric.LinearAlgebra.Matrix.Class where

-- import Numeric.LinearAlgebra.Vector

class Functor m => Matrix m where
  mDim     :: m a -> Int
  mElement :: m a -> Int -> Int -> a
  mIndexOf :: (Ord a) => (a -> a -> Bool) -> m a -> (Int, Int)
  mZip     :: (a -> b -> c) -> m a -> m b -> m c
  -- | mFold is foldl1'
  mFold    :: (a -> a -> a) -> m a -> a
  det      :: Num a => m a -> a

{-# INLINE mApply #-}
mApply :: Functor f => f (a -> b) -> a -> f b
mApply f m = fmap ($ m) f

(.+.) :: (Num k, Matrix m) => m k -> m k -> m k
(.+.) = mZip (+)
(.-.) :: (Num k, Matrix m) => m k -> m k -> m k
(.-.) = mZip (-)

{-
(.*.) :: (k ~ Element m, Num k, Matrix m) => m -> m -> m
m .*. n = mIdxMap m $ \i j -> sum [ mElement m i k * mElement n k j | k <- [ 0 .. 3 ] ]
-}
{-
(.*>) :: (k ~ Element m, k ~ Scalar v, Num k, Matrix m, Vector v)
      => m -> v -> v
m .*> v | vDim v == mDim m = flip vIdxMap v $ \k -> sum [ mElement m i k * vElement v k | i <- [ 0 .. mDim m ] ]
        | otherwise        = error "Dimensions do not match"
-}
-- (*.) :: (k ~ Element m, Num k, Matrix m) => m -> k -> m
-- m *. k = mMap (k*) m

{-
transpose :: Matrix m => m -> m
transpose m = flip mIdxMap m $ \i j -> mElement j i m
-}

max_stars_repo_name: dagit/lin-alg

max_issues_repo_name: dagit/lin-alg

max_forks_repo_name: dagit/lin-alg


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeOperators         #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}

{-|
Module      : Numeric.LinearAlgebra.MatrixFree
Description : Represent matrices as a pair of functions for forward and adjoint.
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module Numeric.LinearAlgebra.MatrixFree
    ( LinearMap(..)
    , (##>)
    , (*#)
    , (+#)
    , (<#>)
    , avgRepeatedFree
    , exactDimsFree
    , eyeFree
    , removeAvg
    , fromL
    , innerDimsFree
    , toL
    , konst
    , diag
    , forward
    , adjoint
    )
where

import           Data.Proxy                     ( Proxy(..) )
import           Data.Type.Equality             ( (:~:)(Refl) )
import           GHC.TypeLits
import           Numeric.LinearAlgebra          ( Transposable(..) )
import           Numeric.LinearAlgebra.Static hiding (konst, diag)
import qualified Numeric.LinearAlgebra.Static  as LS
import qualified Data.Vector.Storable          as V

-- | A linear operator/map between two vector spaces
data LinearMap (m :: Nat) (n :: Nat) where
  LinearMap :: (KnownNat m, KnownNat n) =>
               (R n -> R m) -- ^ Forward Function
            -> (R m -> R n) -- ^ Adjoint FunctionLS
            -> LinearMap m n

-- | Extract the foward function from a `LinearMap`
forward :: LinearMap m n -> (R n -> R m)
forward (LinearMap f _) = f

-- | Extract the adjoint function from a `LinearMap`
adjoint :: LinearMap m n -> (R m -> R n)
adjoint (LinearMap _ a) = a

instance Transposable (LinearMap n m) (LinearMap m n) where
   -- We do not have different functions for conjugate transposes in the
   -- LinearMap constructor (yet?)
    tr  (LinearMap f a) = LinearMap a f
    tr' = tr

---------------------------------------------------------------------
--        Conversions between matrices and matrix free forms       --
---------------------------------------------------------------------

-- | Convert a static matrix into the matrix free data type. This merely
-- wraps the matrix in a function (eta abstraction); it does not determine
-- the function associated with a given matrix.
fromL :: (KnownNat m, KnownNat n) =>
         L m n -- ^ Static matrix to convert to matrix free data type
      -> LinearMap m n
fromL a = LinearMap (a #>) (tr a #>)

-- | Convert a matrix free linear map into a static matrix.
--
-- This methods is very inefficient as it traverses several components
-- linearly several times.
toL :: (KnownNat m, KnownNat n) => LinearMap m n -> L m n
toL (LinearMap f _) =
    matrix $ concatMap (V.toList . extract . f) (toColumns eye)

---------------------------------------------------------------------------
--                           Example LinearMaps                          --
---------------------------------------------------------------------------

-- | The identity matrix in matrix free form.
eyeFree :: (KnownNat n) => LinearMap n n
eyeFree = LinearMap id id

-- | Take the average of a vector.
avg :: forall n. (KnownNat n) => R n -> Double
avg v = v <.> 1 / n' where n' = fromIntegral . natVal $ (Proxy :: Proxy n)

-- | mean(v) * 1, where 1 is the ones vector. of length v.
avgRepeatedFree :: (KnownNat n) => LinearMap n n
avgRepeatedFree = LinearMap f f where f v = LS.konst (avg v)

-- | Remove the average value of a vector from each element of the vector.
removeAvg :: (KnownNat n) => LinearMap n n
removeAvg = eyeFree +# ((-1) *# avgRepeatedFree)

-- | The constant matrix, which sums the vector and sets every element of
-- the resulting vector to that sum.
konst :: (KnownNat m, KnownNat n) => Double -> LinearMap m n
konst k = LinearMap sumAndRepeat sumAndRepeat
  where
    sumAndRepeat :: (KnownNat p, KnownNat q) => R p -> R q
    sumAndRepeat v = LS.konst (v <.> LS.konst k)

-- | Convert a vector into a diagonal matrix.
diag :: (KnownNat n) => R n -> LinearMap n n
diag r = LinearMap f f
  where
    f v = zipWithVector (*) r v

-- TODO: Define fft in matrix free form using FFT and IFFT

---------------------------------------------------------------------
--                Operators for composing LinearMaps               --
---------------------------------------------------------------------

-- | Apply matrix free LinearMap to a vector
(##>) :: LinearMap n m -> R m -> R n
(##>) (LinearMap f _) = f

infixr 8 ##>

-- | Matrix multiply of two `LinearMap`s
(<#>) :: LinearMap m n
      -> LinearMap n p
      -> LinearMap m p
(<#>) (LinearMap f_a a_a) (LinearMap f_b a_b) =
    LinearMap (f_a . f_b) (a_b . a_a)

infixr 8 <#>

-- | Multiply a LinearMap matrix by a constant real number
(*#) :: (KnownNat m, KnownNat n) => Double -> LinearMap m n -> LinearMap m n
(*#) s (LinearMap f a) = LinearMap (\v -> LS.konst s * f v) (\v -> LS.konst s * a v)

infixr 8 *#

-- | Add two `LinearMap`s together if they have the same shape.
(+#) :: (KnownNat m, KnownNat n) =>
        LinearMap m n
     -> LinearMap m n
     -> LinearMap m n
(+#) (LinearMap f_a a_a) (LinearMap f_b a_b) =
    LinearMap (\v -> f_b v + f_a v) (\v -> a_b v + a_a v)

infixr 7 +#

-----------------------------------------------------------------------------
--                    Runtime Proofs on LinearMap shape                    --
-----------------------------------------------------------------------------

-- | Useful for constraining two dependently typed LinearMap matrices to
-- match each other in dimensions when they are unknown at compile-time.
exactDimsFree :: forall m n j k.
                 (KnownNat n, KnownNat m, KnownNat j, KnownNat k) =>
                 LinearMap m n
              -> Maybe (LinearMap j k)
exactDimsFree (LinearMap f b) = do
    Refl <- sameNat (Proxy :: Proxy m) (Proxy :: Proxy j)
    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy k)
    return $ LinearMap f b

-- | Useful for constraining that two dependently typed LinearMap matrices
-- match in their inner dimensions (enabling matrix multiplication) each
-- other in dimensions when they are unknown at compile-time.
innerDimsFree :: forall m n p q. (KnownNat n, KnownNat p) =>
                 LinearMap m n   -- ^ The left LinearMap in a matrix multiply
              -> LinearMap p q   -- ^ The right LinearMap in a matrix multiply
              -> Maybe (n :~: p) -- ^ Proof that the inner dimensions are equal
innerDimsFree (LinearMap _ _) (LinearMap _ _) = do
    Refl <- sameNat (Proxy :: Proxy n) (Proxy :: Proxy p)
    return (Refl)

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
module NumberParser
  ( parseNumber
  , parseFloat
  , parseComplex
  ) where

import Data.Complex
import LispVal
import Numeric
import Text.ParserCombinators.Parsec

parseNumber :: Parser LispVal
parseNumber = parseDecimal1 <|> parseDecimal2 <|> parseHex <|> parseOct <|> parseBin

parseDecimal1 :: Parser LispVal
parseDecimal1 = Number . read <$> many1 digit

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

parseHex :: Parser LispVal
parseHex = do
  try $ string "#x"
  x <- many1 hexDigit
  return $ Number (hex2dig x)

parseOct :: Parser LispVal
parseOct = do
  try $ string "#o"
  x <- many1 octDigit
  return $ Number (oct2dig x)

parseBin :: Parser LispVal
parseBin = do
  try $ string "#b"
  x <- many1 (oneOf "10")
  return $ Number (bin2dig x)

oct2dig x = fst $ head (readOct x)

hex2dig x = fst $ head (readHex x)

bin2dig = bin2dig' 0

bin2dig' digint "" = digint
bin2dig' digint (x:xs) =
  let old =
        2 * digint +
        (if x == '0'
           then 0
           else 1)
   in bin2dig' old xs

toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromIntegral n

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

parseComplex :: Parser LispVal
parseComplex = do
  x <- try parseFloat <|> parseDecimal1 <|> parseDecimal2
  char '+'
  y <- try parseFloat <|> parseDecimal1 <|> parseDecimal2
  char 'i'
  return $ Complex (toDouble x :+ toDouble y)

max_stars_repo_name: anushriadhia/haskellScheme

max_issues_repo_name: anushriadhia/haskellScheme

max_forks_repo_name: anushriadhia/haskellScheme


text:
{-# LANGUAGE FlexibleInstances #-}

import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import Data.Map.Strict (Map, (!), insert, elems, fromList, toList, findWithDefault, size, empty, member, findMin, findMax, singleton, filter)
import qualified Data.Map.Strict as M
import Data.Set (Set)
import qualified Data.Set as S
--import qualified Data.Array as A
import Data.List (find, intercalate, intersperse, permutations, inits, tails, isPrefixOf)
import Data.List.Split (splitOn)
-- import Data.Complex (Complex((:+)), realPart, imagPart) -- define my own complex

import UI.NCurses
import Data.Char (chr, ord)
-- import Data.Complex (Complex((:+)))
import Data.Maybe (fromJust)
type Instructions = Map Integer Integer

data ComputerState = Ready | Running | Blocked | Done
  deriving (Show, Eq)

data Computer = Computer {
  state :: ComputerState,
  memory :: Map Integer Integer,
  iptr :: Integer, -- instruction pointer
  base :: Integer, -- base offset
  input :: [Integer],
  output :: [Integer]}

computer0 = Computer { state = Ready, memory = empty, iptr = 0, base = 0, input = [], output = [] }

run :: Computer -> Computer -- (Integer, Integer) -> Instructions -> [Integer] -> [Integer]
run c
  | state c == Ready = run $ c {state = Running}
  | state c == Done = c
  | state c == Blocked = if input c == [] then c else run $ c {state = Running}
  | otherwise = -- running
    --traceShow (iptr c, base c, memory c) $
    case instr `mod` 100 of
      1 -> -- add
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 + arg 2) instructions}
      2 -> -- multiply
        run $ c {iptr = i+4, memory = insert (addr 3) (arg 1 * arg 2) instructions}
      3 -> -- read input
        if (null $ input c)
        then c {state = Blocked}
        else run $ c {iptr = i+2, memory = insert (addr 1) (head $ input c) instructions, input = tail $ input c}
      4 -> -- output
        --traceShow (arg 1) $
        run $ c {iptr = i+2, output = output c ++ [arg 1]}
      5 -> -- jump-if-true
        run $ c {iptr = if arg 1 == 0 then i+3 else arg 2}
      6 -> -- jump-if-false
        run $ c {iptr = if arg 1 == 0 then arg 2 else i+3}
      7 -> -- less than
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 < arg 2 then 1 else 0) instructions}
      8 -> -- equals
        run $ c {iptr = i+4, memory = insert (addr 3) (if arg 1 == arg 2 then 1 else 0) instructions}
      9 -> -- set relative base
        run $ c {iptr = i+2, base = base c + arg 1}
      99 -> -- halt
        run $ c {state = Done}
      _ -> error "unknown opcode"
  where instructions = memory c
        i = iptr c
        instr = instructions!i
        ii x = findWithDefault 0 x instructions
        arg :: Integer -> Integer
        arg n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ ii $ i+n
          1 -> ii $ i+n
          2 -> ii ((ii $ i+n) + base c)
          _ -> error "bad argument mode"
        addr n = case (instr `mod` (100*10^n)) `div` (10*10^n) of
          0 -> ii $ i+n
          1 -> error "address in mode 1"
          2 -> --trace "address in mode 2" $
            ii (i+n) + base c
          _ -> error $ "address in unknown mode"

instance (Show a) => Show (Complex a) where
  show (C x y) = "C " ++ show x ++ " " ++ show y

instance (Eq a) => Eq (Complex a) where
  (C a ai) == (C b bi) = a == b && ai == bi

instance (Eq a, Ord a) => Ord (Complex a) where
  compare (C a ai) (C b bi) | a < b || a == b && ai < bi = LT
                            | a > b || a == b && ai > bi = GT
                            | otherwise = EQ

-- Note that we are using a plane where imaginary (y-axis) numbers
-- increase downward. I think this means that right and left are
-- reversed.

data Complex a = C a a

-- instance Num (Complex Int) where
--   abs c = error "Can't take abs of Complex Int"

instance Num a => Num (Complex a) where
  (C x y) + (C u v) = C (x+u) (y+v)
  (C x y) * (C u v) = C (x*u-y*v) (x*v+y*u)
  fromInteger n = C (fromInteger n) 0
  abs (C x y) = C (abs x) y -- sqrt $ x*x + y*y
  signum (C x y) = C (signum x) 0
  negate (C x y) = C (negate x) (negate y)

type Point = Complex Int
type Direction = Complex Int
type Scaffold = Set Point
data Move = R | L | F Int | Err
  deriving (Show, Eq)

show' :: [Move] -> String
show' ms = intercalate "," $ map showMove ms
  where showMove R = "R"
        showMove L = "L"
        showMove (F n) = show n
        showMove Err = error "showing err"

charToDirection '>' = C 1 0
charToDirection '<' = C (-1) 0
charToDirection '^' = C 0 (-1)
charToDirection 'v' = C 0 1
charToDirection c = error $ "bad direction" ++ [c]

fromJust' :: String -> Maybe a -> a
fromJust' s Nothing = error s
fromJust' s (Just x) = x

path :: Scaffold -> (Point, Direction) -> [Move]
path scaffold (pt, dir)
  | S.size scaffold < 2 = []
  -- | otherwise = (traceShowId $ rotation ++ [F steps]) ++ path scaffold' (pt', dir')
  | otherwise = rotation ++ [F steps] ++ path scaffold' (pt', dir')
  where left :: Complex Int
        left = C 0 (-1)
        dir' = fromJust' "dir'" $ find (\d -> (pt + d) `S.member` scaffold) [left^i | i <- [1 .. 4]]
        directions = [(left, [L]), (left^2, [L,L]), (left^3, [R]), (left^4, [Err])]
        rotation = snd $ fromJust' ("rotation: "++ show (dir', directions))  $
          find (\(d, m) -> dir * d == dir') directions
        (pt', steps, scaffold') = walk (pt, 0, scaffold)
        walk (p, n, s) | not $ (p+dir') `S.member` s = (p, n, s)
                       | (p+dir'*left) `S.member` s = walk (p+dir', n+1, s) -- intersection, don't delete
                       | otherwise = walk (p+dir', n+1, S.delete p s)

testIntersection scaffold pt = all (`S.member` scaffold) [ pt + d | d <- [l, l*l, l*l*l, l*l*l*l]]
  where l = C 0 1

data Function = A | B | C' deriving (Show)

substitute :: [Move] -> [Function]
-- substitute (L: F n: ms) = A: (substitute $ F (n-8): ms)
-- substitute (F n: ms) = B n: substitute ms
-- substitute (R: F n: ms) = C': (substitute $ F (n-8): ms)
-- substitute [] = []
-- substitute ms = error $ "substitute: can't do " ++ show ms
substitute (L: F 12: L: F 12: R: F 12: ms) = A: substitute ms
substitute (L: F 8: L: F 8: R: F 12: L: F 8: L: F 8: ms) = B: substitute ms
substitute (L: F 10: R: F 8: R: F 12: ms) = C': substitute ms
substitute [] = []
substitute ms = error $ "substitute: can't do " ++ show ms

main = do
  -- [instructionFile] <- getArgs
  instructionStrings <- readFile "17.input.txt" -- instructionFile
  let instructions = fromList . zip [0 ..] $ map read $ splitOn "," instructionStrings

  putStrLn "Part 1"
  let c = run $ computer0{memory = instructions}
  let view = output c
  let viewStr = map (chr . fromIntegral) view
  putStrLn $ viewStr
  let listView = [(C i j, c) | (j, line) <- zip [0 ..] (lines viewStr)
                              , (i, c) <- zip [0 ..] line
                              , c `elem` "#<>^vV"]

  let scaffold = S.fromList $ map fst listView

  let intersections = S.filter (testIntersection scaffold) scaffold
  print $ S.foldl' (\s (C x y) -> s+x*y) 0 intersections

  putStrLn "Part 2"
  let (r, d) = fromJust' "find robot" $ find (\(pt, v) -> v `elem` "<>^vV") listView
  let p = path scaffold (r, charToDirection d)
  -- print p
  -- let cs = substitute p
  -- print cs
  -- print $ length cs
  let ls = map (fromIntegral . ord) $ unlines
        [ "A,A,B,C,C,A,B,C,A,B"
        , "L,12,L,12,R,12"
        , "L,8,L,8,R,12,L,8,L,8"
        , "L,10,R,8,R,12"
        , "n"
        ]
  print ls
  -- let c = run $ computer0 { memory = insert 0 2 instructions
  --                         , input = ls }

  -- print $ output c

  print "Part 2 - alt"
  print p
  print $ length p
  let substrings xs = concatMap inits $ tails xs
  let functions = Prelude.filter (\l -> length l > 1 && (length $ show' l) < 21) $ substrings p
  print $ length $ functions
  -- let functions3 = [(a,b,c) | a <- functions, b <- functions, show' a < show' b, c <- functions, show' b < show' c]
  let f xs ys = show' xs < show' ys
  let functions3 = [[a,b,c] | a <- inits $ take 10 p, length a > 0, c <- functions, b <- Prelude.filter (f c) functions]
  print $ length $ functions3

  -- print $ Prelude.filter (\[a,b,c] -> (length $ show' a) == 14 && (length $ show' b) == 13 && (length $ show' c) == 20) functions3
  let test3 [] fs = True
      test3 pth fs = any id [test3 (drop (length f) pth) fs | f <- fs, isPrefixOf f pth]
  let goodFunctions = Prelude.filter (test3 p) $ functions3
  let goodFunction = head $ goodFunctions
  print $ goodFunction
  let programs [] fs = [[]]
      programs pth fs = [ (l: ls) | (l, f) <- zip ['A' .. 'C'] fs
                                  , f `isPrefixOf` pth
                                  , ls <- programs (drop (length f) pth) fs
                                  ]
  let program = intersperse ',' $ head $ programs p goodFunction
  let functions = map show' goodFunction
  print program
  print functions
  let ls = unlines ((program: functions) ++ ["n"])
  print ls
  print $ map (fromIntegral . ord) ls
  let c = run $ computer0 { memory = insert 0 2 instructions
                          , input = map (fromIntegral . ord) ls }

  print $ last $ output c
  print $ map (chr . fromIntegral) $ init $ output c

-- After first run:
-- [A,B 4,
--  A,B 4,
--  C',B 4,

--  A,B 4,
--  A,B 4,
--  C',B 4,

--  A,B 0,
--  A,B 0,
--  C',B 4,
--  A,B 0,
--  A,B 0,

--  A,B 2,
--  C',B 0,
--  C',B 4,

--  A,B 2,
--  C',B 0,
--  C',B 4,

--  A,B 4,
--  A,B 4,
--  C',B 4,

--  A,B 0,
--  A,B 0,
--  C',B 4,
--  A,B 0,
--  A,B 0,

--  A,B 2,
--  C',B 0,
--  C',B 4,

--  A,B 4,
--  A,B 4,
--  C',B 4,

--  A,B 0,
--  A,B 0,
--  C',B 4,
--  A,B 0,
--  A,B 0]


text:
module Main where

import Math.Probably.EMPCA
import Math.Probably.PCA
import Math.Probably.Datasets

import Data.Random
import System.Environment
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Data
import qualified Data.Vector.Storable as VS

main = do
  fnm:_ <- getArgs
  Right plants <- iris fnm
  let plantvs = map (VS.fromList . plantToVector) plants
      basis = stat $ fromRows $ plantvs
      (enc, _) = pcaN 3 basis
      covar = cov $ fromRows $ map enc plantvs
      covarPre = cov $ fromRows $ plantvs
  disp 2 $ covarPre
  disp 2 $ covar
  em <- sample $ emPca 3 10 (Just cinit) plantvs
  let covarEm = cov $ fromRows $ map (applyEmPca em) plantvs
  disp 2 $ covarEm

cinit = (4><3)
        [ 0.6794,    0.0841,    0.6963,
          0.0731,    0.5255,    0.2313,
          0.6621,    0.0939,    0.7125,
          0.5875,    0.4601,    0.0754  ]

--  mapM_ print plants

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
-----------------------------------------------------------------------------
-- |
-- Module  :  ForSyDe.DFT
-- Copyright   :  (c) ForSyDe Group, KTH 2007-2008
-- License     :  BSD-style (see the file LICENSE)
--
-- Maintainer  :  forsyde-dev@ict.kth.se
-- Stability   :  experimental
-- Portability :  portable
--
-- This module includes the standard Discrete Fourier Transform (DFT)
-- function, and a fast Fourier transform (FFT) algorithm, for
-- computing the DFT, when the input vectors' length is a power of 2.
-----------------------------------------------------------------------------
module ForSyDe.Shallow.Utility.DFT(dft, fft) where

import ForSyDe.Shallow.Core.Vector
import Data.Complex

-- | The function 'dft' performs a standard Discrete Fourier Transformation
dft :: Int -> Vector (Complex Double) -> Vector (Complex Double)
dft bigN x | bigN == (lengthV x) = mapV (bigX_k bigN x) (nVector x)
           | otherwise = error "DFT: Vector has not the right size!"
  where
    nVector x'        = iterateV (lengthV x') (+1) 0
    bigX_k bigN' x' k = sumV (zipWithV (*) x' (bigW' k bigN'))
    bigW' k' bigN'    = mapV (** k') (mapV cis (fullcircle bigN'))
    sumV              = foldlV (+) (0:+ 0)

fullcircle :: Int -> Vector Double
fullcircle n = fullcircle1 0 (fromIntegral n) n
  where
    fullcircle1 l m n'
      | l == m    = NullV
      | otherwise = -2*pi*l/(fromIntegral n')
                    :> fullcircle1 (l+1) m n'

-- | The function 'fft' implements a fast Fourier transform (FFT)
-- algorithm, for computing the DFT, when the size N is a power of 2.
fft :: Int -> Vector (Complex Double) -> Vector (Complex Double)
fft bigN xv | bigN == (lengthV xv) = mapV (bigX xv) (kVector bigN)
            | otherwise = error "FFT: Vector has not the right size!"

kVector :: (Num b, Num a, Eq a) => a -> Vector b
kVector bigN = iterateV bigN (+1) 0

bigX :: Vector (Complex Double) -> Int -> Complex Double
bigX (x0:>x1:>NullV) k | even k = x0 + x1 * bigW 2 0
                       | odd k  = x0 - x1 * bigW 2 0
bigX xv k = bigF_even k + bigF_odd k * bigW bigN (fromIntegral k)
  where bigF_even k' = bigX (evens xv) k'
        bigF_odd k' = bigX (odds xv) k'
        bigN = lengthV xv

bigW :: Int -> Int -> Complex Double
bigW bigN k = cis (-2 * pi * (fromIntegral k) / (fromIntegral bigN))

evens :: Vector a -> Vector a
evens NullV   = NullV
evens (v1:>NullV) = v1 :> NullV
evens (v1:>_:>v)  = v1 :> evens v

odds :: Vector a -> Vector a
odds NullV    = NullV
odds (_:>NullV)   = NullV
odds (_:>v2:>v)   = v2 :> odds v

max_stars_repo_name: hanikat/forsyde-shallow

max_issues_repo_name: hanikat/forsyde-shallow

max_forks_repo_name: hanikat/forsyde-shallow


text:
{-# LANGUAGE Rank2Types #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Packed.Tri
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Triangular views of packed matrices.
--

module Numeric.LinearAlgebra.Packed.Tri (
    -- * Immutable interface

    -- ** Vector multiplication
    triMulVector,

    -- ** Vector solving
    triSolvVector,

    -- * Mutable interface
    triCreate,

    -- ** Vector multiplication
    triMulVectorM_,

    -- ** Vector solving
    triSolvVectorM_,

    ) where

import Control.Monad( when )
import Control.Monad.ST( ST, runST, unsafeIOToST )
import Text.Printf( printf )

import Numeric.LinearAlgebra.Vector( Vector, STVector )
import qualified Numeric.LinearAlgebra.Vector as V
import Numeric.LinearAlgebra.Packed.Base( Packed, RPacked, STPacked )
import qualified Numeric.LinearAlgebra.Packed.Base as P
import Numeric.LinearAlgebra.Types
import qualified Foreign.BLAS as BLAS

-- | A safe way to create and work with a mutable Tri Packed before returning
-- an immutable one for later perusal.
triCreate :: (Storable e)
          => (forall s. ST s (Tri (STPacked s) e))
          -> Tri Packed e
triCreate mt = runST $ do
    (Tri u d ma) <- mt
    a <- P.unsafeFreeze ma
    return $ Tri u d a

-- | @triMulVector trans a x@ returns @op(a) * x@, where @op(a)@ is
-- determined by @trans@.
triMulVector :: (BLAS2 e)
             => Trans
             -> Tri Packed e
             -> Vector e
             -> Vector e
triMulVector trans a x =
    V.create $ do
        x' <- V.newCopy x
        triMulVectorM_ trans a x'
        return x'

-- | @triMulVectorM_ a x@ sets @x := op(a) * x@, where @op(a)@ is determined
-- by @trans@.
triMulVectorM_ :: (RPacked p, BLAS2 e)
               => Trans -> Tri p e
               -> STVector s e
               -> ST s ()
triMulVectorM_ trans (Tri uplo diag a) x = do
    na <- P.getDim a
    nx <- V.getDim x
    let n = nx

    when (nx /= n) $ error $
        printf ("triMulVectorM_"
                ++ " _"
                ++ " (Tri _ _ <packed matrix with dim %d>)"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch")
               na
               nx

    unsafeIOToST $
        P.unsafeWith a $ \pa ->
        V.unsafeWith x $ \px ->
            BLAS.tpmv uplo trans diag n pa px 1

-- | @triSolvVector trans a x@ returns @op(a) \\ x@, where @op(a)@ is
-- determined by @trans@.
triSolvVector :: (BLAS2 e)
             => Trans
             -> Tri Packed e
             -> Vector e
             -> Vector e
triSolvVector trans a x =
    V.create $ do
        x' <- V.newCopy x
        triSolvVectorM_ trans a x'
        return x'

-- | @triSolvVectorM_ a x@ sets @x := op(a) \\ x@, where @op(a)@ is determined
-- by @trans@.
triSolvVectorM_ :: (RPacked p, BLAS2 e)
               => Trans -> Tri p e
               -> STVector s e
               -> ST s ()
triSolvVectorM_ trans (Tri uplo diag a) x = do
    na <- P.getDim a
    nx <- V.getDim x
    let n = nx

    when (nx /= n) $ error $
        printf ("triMulVectorM_"
                ++ " _"
                ++ " (Tri _ _ <packed matrix with dim %d>)"
                ++ " <vector with dim %d>"
                ++ ": dimension mismatch")
               na
               nx

    unsafeIOToST $
        P.unsafeWith a $ \pa ->
        V.unsafeWith x $ \px ->
            BLAS.tpsv uplo trans diag n pa px 1

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}

module HVX.Internal.Primitives
  ( Var
  , Fun(..)
  , Expr(..)
  , getFun
  , getJacobian
  , getProperties
  ) where

import Numeric.LinearAlgebra

import HVX.Internal.Matrix
import HVX.Internal.DCP

type Var = String

data Expr vex mon where
  EConst :: Mat -> Expr 'Affine 'Const
  EVar   :: Var -> Expr 'Affine 'Nondec
  EFun   ::
#ifndef DISABLE_EXPR_CXT
      ValidVex vex =>
#endif
      Fun v1 m1 -> Expr v2 m2 -> Expr vex mon
  EAdd   ::
#ifndef DISABLE_EXPR_CXT
      ValidVex vex =>
#endif
    Expr v1 m1 -> Expr v2 m2 -> Expr vex mon

instance Show (Expr vex mon) where
  show (EConst mat) = show mat
  show (EVar   var) = var
  show (EFun   f a) = show f ++ "(" ++ show a ++ ")"
  show (EAdd   a b) = show a ++ " + " ++ show b

getProperties :: (GetVex vex, GetMon mon) => Expr vex mon -> String
getProperties expr = getVex expr ++ " " ++ getMon expr

data Fun (vex :: Vex) (mon :: Mon) where
  Mul                :: Mat -> Fun 'Affine 'Nonmon
  Abs                :: Fun 'Convex 'Nonmon
  Neg                :: Fun 'Affine 'Noninc
  Log                :: Fun 'Concave 'Nondec
  Exp                :: Fun 'Convex 'Nondec
  LogSumExp          :: Fun 'Convex 'Nondec
  Max                :: Fun 'Convex 'Nondec
  Min                :: Fun 'Concave 'Nondec
  Norm               :: Double -> Fun 'Convex 'Nonmon
  Berhu              :: Double -> Fun 'Convex 'Nonmon
  Huber              :: Double -> Fun 'Convex 'Nonmon
  Quadform           :: Mat -> Fun 'Convex 'Nonmon
  PowBaseP0          :: Fun 'Affine 'Const
  PowBaseP01         :: Double -> Fun 'Concave 'Nondec
  PowBaseP1          :: Fun 'Affine 'Nondec
  PowBaseP1InfEven   :: Integer -> Fun 'Convex 'Nonmon
  PowBaseP1InfNotInt :: Double -> Fun 'Convex 'Nondec

instance Show (Fun vex mon) where
  show (Mul _) = "mul"
  show Abs = "abs"
  show Neg = "-"
  show Log = "log"
  show Exp = "exp"
  show LogSumExp = "log_sum_exp"
  show Max = "max"
  show Min = "min"
  show (Norm p) = "norm" ++ show p
  show (Berhu m) = "berhu" ++ show m
  show (Huber m) = "huber" ++ show m
  show (Quadform m) = "quadform" ++ show m
  show PowBaseP0 = "pow0"
  show (PowBaseP01 p) = "pow" ++ show p
  show PowBaseP1 = "pow1"
  show (PowBaseP1InfEven p) = "pow" ++ show p
  show (PowBaseP1InfNotInt p) = "pow" ++ show p

getFun :: Fun vex mon -> Mat -> Mat
getFun (Mul a) x = a <> x
getFun Abs x = abs x
getFun Neg x = negate x
getFun Exp x = exp x
getFun Log x
  | anyMat (<= 0) x = error "Cannot take log of number <= 0."
  | otherwise = log x
getFun LogSumExp x = reduceMat (log . sum) $ exp x
getFun Max x = reduceMat maximum x
getFun Min x = reduceMat minimum x
getFun (Norm p) x = scalarMat $ lpnorm p x
getFun (Berhu m) x =
  cmap (\y -> if abs y <= m then abs y else (abs y ** 2 + m ** 2) / (2 * m)) x
getFun (Huber m) x =
  cmap (\y -> if abs y <= m then abs y ** 2 else 2 * m * abs y - m ** 2) x
getFun (Quadform m) x = tr x <> m <> x
getFun PowBaseP0 x = konst 1.0 ((rows x), (cols x))
getFun (PowBaseP01 p) x
  | allMat (0 <=) x = matrixPow p x
  | otherwise = error "Cannot raise negative number to power p with 0 < p < 1."
getFun PowBaseP1 x = x
getFun (PowBaseP1InfEven p) x = matrixPow (fromIntegral p) x
getFun (PowBaseP1InfNotInt p) x
  | allMat (0 <=) x = matrixPow p x
  | otherwise = error "Cannot raise negative number to a nonintegral power of p for p > 1."

getJacobian :: Fun vex mon -> Mat -> Mat
getJacobian (Mul a) _ = a
getJacobian Abs x = diagMat $ signum x
getJacobian Neg x = (-1) * ident (rows x)
getJacobian Log x = diagMat $ 1 / x
getJacobian Exp x = diagMat $ exp x
getJacobian LogSumExp x = tr $ scale (1 / sumElements (exp x)) $ exp x
getJacobian Max x = tr $ ei (rows x) i
  where (i, _) = maxIndex x
getJacobian Min x = tr $ ei (rows x) i
  where (i, _) = minIndex x
getJacobian (Norm p) x = tr $ if x == zeroVec n
  then
    zeroVec n
  else
    scale (1 / denominator) numerator
  where denominator = lpnorm p x ** (p - 1)
        numerator = x * abs x ** (pMat - 2)
        pMat = (1><1) [p]
        n = rows x
getJacobian (Berhu m) x = diagMat $
  cmap (\y -> if abs y <= m then signum y else y / m) x
getJacobian (Huber m) x = diagMat $
  cmap (\y -> if abs y <= m then 2 * y else 2 * m * signum y) x
getJacobian (Quadform m) x = scale 2 $ tr x <> m
getJacobian PowBaseP0 x = konst 0.0 (n, n)
  where n = rows x
getJacobian (PowBaseP01 p) x = scale p $ diagMat $ matrixPow (p - 1) x
getJacobian PowBaseP1 x = ident n
  where n = rows x
getJacobian (PowBaseP1InfEven intP) x = scale p $ diagMat $ matrixPow (p - 1) x
  where p = fromIntegral intP
getJacobian (PowBaseP1InfNotInt p) x = scale p $ diagMat $ matrixPow (p - 1) x

max_stars_repo_name: chrisnc/hvx

max_issues_repo_name: chrisnc/hvx

max_forks_repo_name: chrisnc/hvx


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE BangPatterns #-}
module Taiji.Utils
    ( BBIndex
    , openBBs
    , queryBB
    , readExpression
    , lp

    , readGenesValidated

    , computeRAS
    , computeSS
    , computeCDF

    , module Taiji.Utils.Matrix
    , module Taiji.Utils.Clustering
    ) where

import Data.BBI.BigBed
import Bio.Data.Experiment
import Bio.Data.Bed
import Bio.Utils.Functions (scale)
import           Data.CaseInsensitive    (mk)
import qualified Data.HashMap.Strict as HM
import qualified Data.Matrix.Unboxed as MU
import qualified Data.Matrix as Mat
import qualified Data.ByteString.Char8 as B
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector.Unboxed.Mutable as UM
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as VM
import System.Random.MWC.Distributions
import System.Random.MWC
import Statistics.Sample (mean, varianceUnbiased)
import Bio.RealWorld.GENCODE (readGenes, Gene(..))

import Taiji.Utils.Clustering
import Taiji.Utils.Matrix
import qualified Taiji.Utils.DataFrame as DF
import Taiji.Prelude

readGenesValidated :: FilePath -> IO [Gene]
readGenesValidated = fmap validate . readGenes
  where
    validate xs
        | null xs = error "No gene was found in the annotation file"
        | all (null . geneTranscripts) xs = error "No transcript was found in the annotation file"
        | otherwise = xs
{-# INLINE readGenesValidated #-}

-- | Open bigbed files.
openBBs :: [(B.ByteString, Maybe (File '[] 'BigBed))]
        -> IO BBIndex
openBBs xs = fmap (HM.fromList . catMaybes) $ forM xs $ \(chr, x) -> case x of
    Nothing -> return Nothing
    Just fl -> do
        bb <- openBBedFile $ fl^.location
        return $ Just (chr, bb)

queryBB :: BEDLike b => b -> BBIndex -> ConduitT () TFBS IO ()
queryBB bed idx = case HM.lookup (bed^.chrom) idx of
    Nothing -> return ()
    Just i -> query (bed^.chrom, bed^.chromStart, bed^.chromEnd) i .| mapC f
  where
    f (chr, s, e, rest) = BEDExt (asBed chr s e) site
      where
        site = SiteInfo
            { _tf_name = mk $ head $ B.split '+' f1
            , _site_affinity = toSiteAffinity $ readInt f2
            , _peak_affinity = toPeakAffinity 100 }
        (f1:f2:_) = B.split '\t' rest

lp :: Int -> [Double] -> Double
lp p = (**(1/fromIntegral p)) . foldl' (+) 0 . map (**fromIntegral p)
{-# INLINE lp #-}

-- | Read RNA expression data
readExpression :: Double    -- ^ Threshold to call a gene as non-expressed
               -> B.ByteString  -- ^ cell type
               -> FilePath
               -> IO (HM.HashMap GeneName (Double, Double)) -- ^ absolute value and z-score
readExpression thres ct fl = do
    c <- B.readFile fl
    let ((_:header):dat) = map (B.split '\t') $ B.lines c
        rowNames = map (mk . head) dat
        dataTable = map (U.fromList . map ((+pseudoCount) . readDouble) . tail) dat
        dataTable' = MU.fromRows $ zipWith U.zip dataTable $
            map computeZscore dataTable :: MU.Matrix (Double, Double)
        idx = fromMaybe (error $ "Cell type:" ++ B.unpack ct ++ " not found!") $
            lookup ct $ zip header [0..]
    return $ HM.fromList $ zip rowNames $ U.toList $ dataTable' `MU.takeColumn` idx
  where
    pseudoCount = 0.1
    computeZscore xs
        | U.length xs == 1 = U.map log xs
        | U.all (<thres) xs = U.replicate (U.length xs) 0
        | U.length xs == 2 = let fc = log $ U.head xs / U.last xs
                             in U.fromList [fc, negate fc]
        | U.all (== U.head xs) xs = U.replicate (U.length xs) 0
        | otherwise = scale xs
{-# INLINE readExpression #-}

computeRAS :: FilePath -> IO (U.Vector Double)
computeRAS fl = do
    mat <- mkSpMatrix readInt fl
    fmap accScore $ runResourceT $ runConduit $
        streamRows mat .| sink (_num_col mat)
  where
    sink n = do
        v <- lift $ UM.replicate n 0.1
        mapM_C $ \(_, xs) -> forM_ xs $ \(i, x) ->
            UM.unsafeModify v (+fromIntegral x) i
        lift $ U.unsafeFreeze v
    accScore xs = U.map (\x -> x * 1000000 / s) xs
      where
        s = U.sum xs

-- | Compute Specificity Score (SS).
computeSS :: DF.DataFrame Double -> DF.DataFrame Double
computeSS df = df{ DF._dataframe_data = mat}
  where
    mat = Mat.fromRows $ map (q_t . normalize) $ Mat.toRows $
        DF._dataframe_data df
    q_t xs = V.map (\x -> e - logBase 2 x) xs
      where
        e = negate $ V.sum $ V.map (\p -> p * logBase 2 p) xs
    normalize xs = V.map (/s) xs
      where
        s = V.sum xs
{-# INLINE computeSS #-}

computeCDF :: DF.DataFrame Double -> IO (V.Vector Double, Double, Double)
computeCDF df = do
    gen <- create
    let std = let rows = filter ((>1) . V.maximum) $ Mat.toRows $ DF._dataframe_data df
                  f xs = (xs, 1 / entropy (normalize xs))
                  n = truncate $ (0.3 :: Double) * fromIntegral (length rows)
                  getFold xs = let m = mean xs in V.map (logBase 2 . (/m)) xs
                  entropy xs = negate $ V.sum $ V.map (\p -> p * logBase 2 p) xs
              in sqrt $ varianceUnbiased $ V.concat $ map (getFold . fst) $
                  take n $ sortBy (comparing snd) $ map f $ rows
    runConduit $ replicateMC nSample (V.toList <$> mkSample gen std) .|
        concatC .| mkCDF 0.001
  where
    ncol = Mat.cols $ DF._dataframe_data df
    nSample = 500000
    mkCDF res = do
        vec <- VM.replicate n 0
        mapM_C $ \x -> do
          let i = min (n-1) $ truncate $ x / res
          VM.unsafeModify vec (+1) i
        v <- V.scanl1' (+) <$> V.unsafeFreeze vec
        return (V.map (/(V.last v)) v, res, fromIntegral $ nSample*ncol)
      where
        n = truncate $ 2 * logBase 2 (fromIntegral ncol) / res
    -- | https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2396180/
    normalize xs = V.map (/s) xs
      where
        s = V.sum xs
    mkSample gen !std = do
        folds <- replicateM ncol $ normal 0 std gen
        return $ q_t $ normalize $ V.fromList $ map (\x -> 2**x) folds
    q_t xs = V.map (\x -> e - logBase 2 x) xs
      where
        e = negate $ V.sum $ V.map (\p -> p * logBase 2 p) xs
{-# INLINE computeCDF #-}

max_stars_repo_name: Taiji-pipeline/Taiji-utils

max_issues_repo_name: Taiji-pipeline/Taiji-utils

max_forks_repo_name: Taiji-pipeline/Taiji-utils


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

import qualified Control.Foldl                as F
import           Control.Lens
import qualified Data.Attoparsec.Text         as A
import qualified Data.ByteString.Lazy         as BS
import           Data.Maybe                   (fromJust)
import           Data.Semigroup               ((<>))
import           Data.Serialize
import           Data.Singletons              (SingI)
import           Data.Singletons.Prelude.List (Head, Last)
import           GHC.Exts
import           GHC.TypeLits
import           Grenade
import           Grenade.Recurrent
import qualified Numeric.LinearAlgebra.Static as SA
import           Options.Applicative
import           Pipes
import qualified Pipes.ByteString             as PBS
import           Pipes.Group
import qualified Pipes.Prelude                as P
import qualified Pipes.Prelude.Text           as PT
import           System.IO                    (BufferMode (..), IOMode (..),
                                               hSetBuffering, stdout, withFile)

-- The defininition for our simple recurrent network.
-- This file just trains a network to generate a repeating sequence
-- of 0 0 1.
--
-- The F and R types are Tagging types to ensure that the runner and
-- creation function know how to treat the layers.
type R = Recurrent
type F = FeedForward
type TRNNIn = 4
type TRNNOut = 50
type TFFNOut = 1
type TRNN = R (LSTM TRNNIn TRNNOut)
type TFFN = F (FullyConnected TRNNOut TFFNOut)
-- type NJetInputs = 4
-- type NJetOutputs = 14
-- type JetLayer = Concat ('D1 4) Trivial ('D1 4) TracksLayer

type TracksNetShape = '[ TRNN, TFFN ]

type TracksNet =
  RecurrentNetwork
    TracksNetShape
    '[ 'D1 TRNNIn, 'D1 TRNNOut, 'D1 TFFNOut ]

type TracksInput = RecurrentInputs TracksNetShape

data Opts = Opts Int (Maybe String) LearningParameters

feedForward' :: Parser Opts
feedForward' =
  Opts
    <$> option auto (long "examples" <> short 'e' <> value 40000)
    <*> optional (strOption (long "input" <> short 'i'))
    <*> ( LearningParameters
          <$> option auto (long "train_rate" <> short 'r' <> value 0.01)
          <*> option auto (long "momentum" <> value 0.9)
          <*> option auto (long "l2" <> value 0.0005)
        )

parseJet :: A.Parser [(S ('D1 TRNNIn), Maybe (S ('D1 1)))]
parseJet = do
  nB <- A.double
  tracks <- A.many1 . A.count 4 $ A.char '\t' *> A.double

  return $
    appLast Nothing (Just . S1D $ SA.konst nB)
    $ S1D . SA.fromList <$> tracks

  where
    appLast _ _ []     = []
    appLast _ y [z]    = [(z, y)]
    appLast x y (z:zs) = (z, x) : appLast x y zs

-- can't export this....
trainRecurrentF
  :: (Num (RecurrentInputs layers), SingI (Last shapes), Monad m)
  => LearningParameters
  -> m (RecurrentNetwork layers shapes, RecurrentInputs layers)
  -> F.FoldM
     m
     [(S (Head shapes), Maybe (S (Last shapes)))]
     (RecurrentNetwork layers shapes, RecurrentInputs layers)
trainRecurrentF lps start = F.FoldM step start done
  where
    step (!rnet, !rins) samp =
      return $ trainRecurrent lps rnet rins samp
    done = return

runRecurrentP
  :: Monad m
  => RecurrentNetwork layers shapes
  -> RecurrentInputs layers
  -> Pipe [S (Head shapes)] [S (Last shapes)] m b
runRecurrentP net inps = do
  xs <- await
  let (_, _, ys) = runRecurrentNetwork net inps xs
  yield ys
  runRecurrentP net inps

main :: IO ()
main = do
    Opts nex maybeIn rate <- execParser (info (feedForward' <**> helper) idm)
    hSetBuffering stdout LineBuffering
    start <-
      case maybeIn of
        Nothing -> randomRecurrent
        Just s ->
          withFile s ReadMode $ \h -> do
            res <- decodeLazy <$> BS.hGetContents h
            case res of
              Left err -> error err
              Right x  -> return x

    putStrLn "Training network..."

    let p =
          PT.stdinLn
          >-> P.map (A.parseOnly parseJet)
          >-> P.concat

    (net :: TracksNet, inps :: TracksInput) <-
      F.impurely P.foldM (trainRecurrentF rate $ return start)
      . over (chunksOf 1000) (maps $ \x -> liftIO (putStrLn "trained 1000 more") >> x)
      $ p >-> P.take nex

    withFile "net.out" WriteMode $ \h ->
      runEffect $ PBS.fromLazy (encodeLazy (net, inps)) >-> PBS.toHandle h

    putStrLn ""

    runEffect $
      p
      >-> P.mapM (
        \x -> do
          putStr . show . fromS1D . fromJust . snd $ last x
          putStr ", "
          return x
        )
      >-> P.map (fmap fst)
      >-> runRecurrentP net inps
      >-> P.mapM_ (print . fromS1D . last)

  where
    fromS1D :: KnownNat n => S ('D1 n) -> Double
    fromS1D (S1D r) = head . toList $ SA.unwrap r

max_stars_repo_name: cspollard/bdiscnn

max_issues_repo_name: cspollard/bdiscnn

max_forks_repo_name: cspollard/bdiscnn


text:
module HVX.Internal.TestUtil
  ( fpequalsApprox
  , pairsOf
  , positiveVectors
  , smallEvenPowersGreaterThanOne
  , smallNonIntPowersGreaterThanOne
  , smallVectors
  , vectors
  ) where

import Control.Monad
import Numeric.LinearAlgebra
import Test.QuickCheck hiding ( (><) )

import HVX.Internal.Matrix
import HVX.Internal.Util

fpequalsApprox :: Double -> Double -> Bool
fpequalsApprox a b
  | a == 0 && b == 0           = True
  | max (abs a) (abs b) < 1e-6 = True
  | otherwise                  = abs ( (a - b) / denom) < tolerance
  where
    denom = max (abs a) (abs b)
    tolerance = 5e-2  -- 5% tolerance is why we call it Approx.

pairsOf :: Gen a -> Gen b -> Gen (a, b)
pairsOf a b = do
  x <- a
  y <- b
  return (x, y)

positiveVectors :: Gen Mat
positiveVectors = do
  nonEmptyL <- arbitrary :: Gen (NonEmptyList Double)
  let l = getNonEmpty nonEmptyL
  return $ abs $ (length l><1) l

smallEvenPowersGreaterThanOne :: Gen Double
smallEvenPowersGreaterThanOne =
  liftM fromIntegral $
    (arbitrary :: Gen Integer) `suchThat` \x -> x > 1 && x <= 20 && even x

smallNonIntPowersGreaterThanOne :: Gen Double
smallNonIntPowersGreaterThanOne =
  arbitrary `suchThat` \x -> x > 1 && x <= 20 && not (isInteger x)

smallVectors :: Gen Mat
smallVectors = do
  len <- choose (1, 100)
  l <- sequence [choose (-10, 10) | _ <- [1..len]]
  return $ (len><1) l

vectors :: Gen Mat
vectors = do
  nonEmptyL <- arbitrary :: Gen (NonEmptyList Double)
  let l = getNonEmpty nonEmptyL
  return $ (length l><1) l

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx


text:
{-# LANGUAGE BangPatterns  #-}
{-# LANGUAGE DeriveGeneric #-}

module GraphBuilder (
    shortDistMat,
    shortestDists,
    joinContext,
    EContext(..),
    EShortDists(..),
    theContextsAndDists,
    ) where

import           Control.DeepSeq                   (NFData (..))
import           Control.Monad                     (ap, join, liftM2)
import           Control.Monad.ST
import           Control.Parallel.Strategies
import           Data.Graph.Inductive              hiding (ap)
import qualified Data.Graph.Inductive.PatriciaTree as GP
import           Data.List                         (sortOn)
import qualified Data.Strict.Tuple                 as T hiding (uncurry)
import           Foreign.C                         (CInt)
import           GHC.Generics
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Devel

data EContext =
       EContext
         { vertex   :: {-# UNPACK #-} !Int
         , dstToCtr :: {-# UNPACK #-} !Double
         , dists    :: {-# UNPACK #-} !(Vector Double)
         , ns       :: {-# UNPACK #-} !(Vector Int)
         }
  deriving (Show, Generic)

instance NFData EContext

data EShortDists =
       EShortDists
         { vertex' :: {-# UNPACK #-} !Int
         , dists'  :: {-# UNPACK #-} !(Vector Double)
         }
  deriving (Show, Generic)

instance NFData EShortDists

shortDistMat :: Int -> [EShortDists] -> Matrix Double
shortDistMat size' !sdists =
  let mat = konst 0 (length', size') :: Matrix Double
      length' = length sdists
  in runST $ do
    m <- unsafeThawMatrix mat
    mapM_ (`setRow` m) $ zip sdists [0 ..]
    unsafeFreezeMatrix m
  where
    setRow (e, pt) mat = setMatrix mat pt 0 . asRow . dists' $ e

joinContext :: EContext -> GP.Gr Double Double -> GP.Gr Double Double
joinContext (EContext pt dstCtr edists vs) graph =
  let zipped = zip (toList edists) (toList vs)
  in (zipped, pt, dstCtr, zipped) & graph

matD2 :: Matrix Double -> Matrix Double
matD2 = join pairwiseD2

theContextsAndDists :: (Int, Int) -- ^ Graph neighborhood size and PCA neighborhood size
                    -> Int
                    -> Matrix Double
                    -> ([EContext], Matrix Double, [Int])
theContextsAndDists (n1, n2) bpt mat =
  (makeContexts . theNearest n1 $ mat', extract' . idxs $ idxs'', idxs''')
  where
    makeContexts = zipWith makeContext [0 :: Int ..] . map unzip
    makeContext a (b, c) = EContext a (mat' `atIndex` (a, bpt)) (fromList b) (fromList c)
    mat' = matD2 mat
    idxs'' = map (toEnum . snd) idxs'
    idxs' = (!! bpt) . theNearest n2 $ mat'
    idxs''' = map snd idxs'
    extract' a = mat ?? (Pos a, All)

theNearest :: Int -> Matrix Double -> [[(Double, Int)]]
theNearest n = filterKNearest . findNearestWithInd
  where
    findNearestWithInd = helper1 . helper2
    filterKNearest = map (take n . filter ((/= 0.0) . fst) . uncurry zip)
    helper1 in' = uncurry zip (map toList (T.fst in'), map (map fromEnum . toList) (T.snd in'))
    helper2 = strUnzip . map (liftM2 (T.:!:) sortVector sortIndex) . toRows

shortestDists :: [Int] -> GP.Gr Double Double -> [EShortDists]
shortestDists inds g = parMap rseq (clean `ap` (sortIndDumpDist . makeList)) inds
  where
    clean x' a = EShortDists x' $ fromList a
    sortIndDumpDist = map T.snd . sortOn T.fst
    makeList = map (Prelude.uncurry (T.:!:)) . join . map (take 1 . unLPath) . flip spTree g

-- UTILITY FUNCTION(S) --

strUnzip :: [T.Pair (Vector Double) (Vector CInt)] -> T.Pair [Vector Double] [Vector CInt]
strUnzip = foldr (\((T.:!:) a b) acc -> ((T.:!:) (a : T.fst acc) (b : T.snd acc))) ((T.:!:) [] [])

-- -- For debugging!
-- remember to import Control.Arrow ((&&&)) when using this function.
-- checkSymRowsCols :: (Container Vector a, Eq a, S.Storable a)
--                     => Matrix a
--                     -> Bool
-- checkSymRowsCols = uncurry (==) . (map sumElements . toColumns &&& map sumElements . toRows)

max_stars_repo_name: emmanueldenloye/manifoldRNC

max_issues_repo_name: emmanueldenloye/manifoldRNC

max_forks_repo_name: emmanueldenloye/manifoldRNC


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

module Backprop.Learn.Run (
    consecutives
  , consecutivesN
  , leadings
  , conduitLearn, conduitLearnStoch
  -- * Encoding and decoding for learning
  , oneHot', oneHot, oneHotR
  , SVG.maxIndex, maxIndexR
  ) where

import           Backprop.Learn.Model
import           Control.Monad
import           Control.Monad.Primitive
import           Control.Monad.Trans.Class
import           Control.Monad.Trans.Maybe
import           Data.Bool
import           Data.Conduit
import           Data.Finite
import           Data.Foldable
import           Data.Proxy
import           GHC.TypeNats
import           Numeric.LinearAlgebra.Static
import           Numeric.LinearAlgebra.Static.Vector
import qualified Data.Conduit.Combinators            as C
import qualified Data.Sequence                       as Seq
import qualified Data.Vector.Generic                 as VG
import qualified Data.Vector.Generic.Sized           as SVG
import qualified System.Random.MWC                   as MWC

consecutives :: Monad m => ConduitT i (i, i) m ()
consecutives = void . runMaybeT $ do
    x <- MaybeT await
    go x
  where
    go x = do
      y <- MaybeT await
      lift $ yield (x, y)
      go y

consecutivesN
    :: forall v n i m. (KnownNat n, VG.Vector v i, Monad m)
    => ConduitT i (SVG.Vector v n i, SVG.Vector v n i) m ()
consecutivesN = conseq (fromIntegral n) .| C.concatMap process
  where
    n = natVal (Proxy @n)
    process (xs, ys, _) = (,) <$> SVG.fromList (toList xs)
                              <*> SVG.fromList (toList ys)

leadings
    :: forall v n i m. (KnownNat n, VG.Vector v i, Monad m)
    => ConduitT i (SVG.Vector v n i, i) m ()
leadings = conseq (fromIntegral n) .| C.concatMap process
  where
    n = natVal (Proxy @n)
    process (xs, _, y) = (, y) <$> SVG.fromList (toList xs)

conseq
    :: forall i m. Monad m
    => Int
    -> ConduitT i (Seq.Seq i, Seq.Seq i, i) m ()
conseq n = void . runMaybeT $ do
    xs <- Seq.replicateM n $ MaybeT await
    go xs
  where
    go xs = do
      _ Seq.:<| xs' <- pure xs
      y <- MaybeT await
      let ys = xs' Seq.:|> y
      lift $ yield (xs, ys, y)
      go ys

conduitLearn
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), Monad m)
    => l
    -> LParam_ I l
    -> LState_ I l
    -> ConduitT a b m (LState_ I l)
conduitLearn l p = go
  where
    go s = do
      mx <- await
      case mx of
        Nothing -> return s
        Just x  -> do
          let (y, s') = runLearn_ l p x s
          yield y
          go s'

conduitLearnStoch
    :: (Learn a b l, Backprop b, MaybeC Backprop (LStateMaybe l), PrimMonad m)
    => l
    -> MWC.Gen (PrimState m)
    -> LParam_ I l
    -> LState_ I l
    -> ConduitT a b m (LState_ I l)
conduitLearnStoch l g p = go
  where
    go s = do
      mx <- await
      case mx of
        Nothing -> return s
        Just x  -> do
          (y, s') <- lift $ runLearnStoch_ l g p x s
          yield y
          go s'

-- | What module should this be in?
oneHot'
    :: (VG.Vector v a, KnownNat n)
    => a                -- ^ not hot
    -> a                -- ^ hot
    -> Finite n
    -> SVG.Vector v n a
oneHot' nothot hot i = SVG.generate (bool nothot hot . (== i))

oneHot
    :: (VG.Vector v a, KnownNat n, Num a)
    => Finite n
    -> SVG.Vector v n a
oneHot = oneHot' 0 1

oneHotR :: KnownNat n => Finite n -> R n
oneHotR = vecR . oneHot

-- | Could be in /hmatrix/.
maxIndexR :: KnownNat n => R (n + 1) -> Finite (n + 1)
maxIndexR = SVG.maxIndex . rVec

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn


text:
{-# LANGUAGE FlexibleInstances, ScopedTypeVariables,
    ViewPatterns #-}
module Tests.Distribution (tests) where

import Control.Applicative ((<$), (<$>), (<*>))
import qualified Control.Exception as E
import Data.List (find)
import Data.Typeable (Typeable)
import Numeric.MathFunctions.Constants (m_tiny,m_huge,m_epsilon)
import Numeric.MathFunctions.Comparison
import Statistics.Distribution
import Statistics.Distribution.Beta           (BetaDistribution)
import Statistics.Distribution.Binomial       (BinomialDistribution)
import Statistics.Distribution.CauchyLorentz
import Statistics.Distribution.ChiSquared     (ChiSquared)
import Statistics.Distribution.Exponential    (ExponentialDistribution)
import Statistics.Distribution.FDistribution  (FDistribution,fDistribution)
import Statistics.Distribution.Gamma          (GammaDistribution,gammaDistr)
import Statistics.Distribution.Geometric
import Statistics.Distribution.Hypergeometric
import Statistics.Distribution.Laplace        (LaplaceDistribution)
import Statistics.Distribution.Lognormal      (LognormalDistribution)
import Statistics.Distribution.Normal         (NormalDistribution)
import Statistics.Distribution.Poisson        (PoissonDistribution)
import Statistics.Distribution.StudentT
import Statistics.Distribution.Transform      (LinearTransform)
import Statistics.Distribution.Uniform        (UniformDistribution)
import Statistics.Distribution.Weibull        (WeibullDistribution)
import Statistics.Distribution.DiscreteUniform (DiscreteUniform)
import Test.Tasty                 (TestTree, testGroup)
import Test.Tasty.QuickCheck      (testProperty)
import Test.Tasty.ExpectedFailure (ignoreTest)
import Test.QuickCheck as QC
import Test.QuickCheck.Monadic as QC
import Text.Printf (printf)

import Tests.ApproxEq  (ApproxEq(..))
import Tests.Helpers   (T(..), Double01(..), testAssertion, typeName)
import Tests.Helpers   (monotonicallyIncreasesIEEE,isDenorm)
import Tests.Orphanage ()

-- | Tests for all distributions
tests :: TestTree
tests = testGroup "Tests for all distributions"
  [ contDistrTests (T :: T BetaDistribution        )
  , contDistrTests (T :: T CauchyDistribution      )
  , contDistrTests (T :: T ChiSquared              )
  , contDistrTests (T :: T ExponentialDistribution )
  , contDistrTests (T :: T GammaDistribution       )
  , contDistrTests (T :: T LaplaceDistribution     )
  , contDistrTests (T :: T LognormalDistribution   )
  , contDistrTests (T :: T NormalDistribution      )
  , contDistrTests (T :: T UniformDistribution     )
  , contDistrTests (T :: T WeibullDistribution     )
  , contDistrTests (T :: T StudentT                )
  , contDistrTests (T :: T (LinearTransform NormalDistribution))
  , contDistrTests (T :: T FDistribution           )

  , discreteDistrTests (T :: T BinomialDistribution       )
  , discreteDistrTests (T :: T GeometricDistribution      )
  , discreteDistrTests (T :: T GeometricDistribution0     )
  , discreteDistrTests (T :: T HypergeometricDistribution )
  , discreteDistrTests (T :: T PoissonDistribution        )
  , discreteDistrTests (T :: T DiscreteUniform            )

  , unitTests
  ]

----------------------------------------------------------------
-- Tests
----------------------------------------------------------------

-- Tests for continuous distribution
contDistrTests :: (Param d, ContDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> TestTree
contDistrTests t = testGroup ("Tests for: " ++ typeName t) $
  cdfTests t ++
  [ testProperty "PDF sanity"              $ pdfSanityCheck     t
  ] ++
  [ (if quantileIsInvCDF_enabled t then id else ignoreTest)
  $ testProperty "Quantile is CDF inverse" $ quantileIsInvCDF t
  , testProperty "quantile fails p<0||p>1" $ quantileShouldFail t
  , testProperty "log density check"       $ logDensityCheck    t
  , testProperty "complQuantile"           $ complQuantileCheck t
  ]

-- Tests for discrete distribution
discreteDistrTests :: (Param d, DiscreteDistr d, QC.Arbitrary d, Typeable d, Show d) => T d -> TestTree
discreteDistrTests t = testGroup ("Tests for: " ++ typeName t) $
  cdfTests t ++
  [ testProperty "Prob. sanity"         $ probSanityCheck       t
  , testProperty "CDF is sum of prob."  $ discreteCDFcorrect    t
  , testProperty "Discrete CDF is OK"   $ cdfDiscreteIsCorrect  t
  , testProperty "log probabilty check" $ logProbabilityCheck   t
  ]

-- Tests for distributions which have CDF
cdfTests :: (Param d, Distribution d, QC.Arbitrary d, Show d) => T d -> [TestTree]
cdfTests t =
  [ testProperty "C.D.F. sanity"        $ cdfSanityCheck         t
  , testProperty "CDF limit at +inf"    $ cdfLimitAtPosInfinity  t
  , testProperty "CDF limit at -inf"    $ cdfLimitAtNegInfinity  t
  , testProperty "CDF at +inf = 1"      $ cdfAtPosInfinity       t
  , testProperty "CDF at -inf = 1"      $ cdfAtNegInfinity       t
  , testProperty "CDF is nondecreasing" $ cdfIsNondecreasing     t
  , testProperty "1-CDF is correct"     $ cdfComplementIsCorrect t
  ]

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

-- CDF is in [0,1] range
cdfSanityCheck :: (Distribution d) => T d -> d -> Double -> Bool
cdfSanityCheck _ d x = c >= 0 && c <= 1
  where c = cumulative d x

-- CDF never decreases
cdfIsNondecreasing :: (Distribution d) => T d -> d -> Double -> Double -> Bool
cdfIsNondecreasing _ d = monotonicallyIncreasesIEEE $ cumulative d

-- cumulative d +∞ = 1
cdfAtPosInfinity :: (Distribution d) => T d -> d -> Bool
cdfAtPosInfinity _ d
  = cumulative d (1/0) == 1

-- cumulative d - ∞ = 0
cdfAtNegInfinity :: (Distribution d) => T d -> d -> Bool
cdfAtNegInfinity _ d
  = cumulative d (-1/0) == 0

-- CDF limit at +∞ is 1
cdfLimitAtPosInfinity :: (Param d, Distribution d) => T d -> d -> Bool
cdfLimitAtPosInfinity _ d
  = Just 1.0 == find (>=1) probs
  where
    probs = map (cumulative d)
          $ takeWhile (< (m_huge/2))
          $ iterate (*1.4) 1

-- CDF limit at -∞ is 0
cdfLimitAtNegInfinity :: (Param d, Distribution d) => T d -> d -> Bool
cdfLimitAtNegInfinity _ d
  = Just 0 == find (<=0) probs
  where
    probs = map (cumulative d)
          $ takeWhile (> (-m_huge/2))
          $ iterate (*1.4) (-1)

-- CDF's complement is implemented correctly
cdfComplementIsCorrect :: (Distribution d, Param d) => T d -> d -> Double -> Bool
cdfComplementIsCorrect _ d x
  = 1 - (cumulative d x + complCumulative d x) <= tol
  where
    tol = prec_complementCDF d

-- CDF for discrete distribution uses <= for comparison
cdfDiscreteIsCorrect :: (Param d, DiscreteDistr d) => T d -> d -> Property
cdfDiscreteIsCorrect _ d
  = counterexample (unlines badN)
  $ null badN
  where
    -- We are checking that:
    --
    -- > CDF(i) - CDF(i-e) = P(i)
    --
    -- Apporixmate equality is tricky here. Scale is set by maximum
    -- value of CDF and probability. Case when all proabilities are
    -- zero should be trated specially.
    badN = [ printf "N=%3i    p[i]=%g\tp[i+1]=%g\tdP=%g\trelerr=%g" i p p1 dp ((p1-p-dp) / max p1 dp)
           | i <- [0 .. 100]
           , let p      = cumulative d $ fromIntegral i - 1e-6
                 p1     = cumulative d $ fromIntegral i
                 dp     = probability d i
                 relerr = ((p1 - p) - dp) / max p1 dp
           , p  > m_tiny || p == 0
           , p1 > m_tiny
           , dp > m_tiny
           , relerr > tol
           ]
    tol = prec_discreteCDF d

logDensityCheck :: (ContDistr d) => T d -> d -> Double -> Property
logDensityCheck _ d x
  = not (isDenorm x)
  ==> ( counterexample (printf "density    = %g" p)
      $ counterexample (printf "logDensity = %g" logP)
      $ counterexample (printf "log p      = %g" (log p))
      $ counterexample (printf "eps        = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
      $ or [ p == 0      && logP == (-1/0)
           , p <= m_tiny && logP < log m_tiny
             -- To avoid problems with roundtripping error in case
             -- when density is computed as exponent of logDensity we
             -- accept either inequality
           ,  (ulpDistance (log p) logP <= 32)
           || (ulpDistance p (exp logP) <= 32)
           ])
  where
    p    = density d x
    logP = logDensity d x

-- PDF is positive
pdfSanityCheck :: (ContDistr d) => T d -> d -> Double -> Bool
pdfSanityCheck _ d x = p >= 0
  where p = density d x

complQuantileCheck :: (ContDistr d) => T d -> d -> Double01 -> Property
complQuantileCheck _ d (Double01 p)
  = counterexample (printf "x0 = %g" x0)
  $ counterexample (printf "x1 = %g" x1)
  $ counterexample (printf "abs err = %g" $ abs (x1 - x0))
  $ counterexample (printf "rel err = %g" $ relativeError x1 x0)
  -- We avoid extreme tails of distributions
  --
  -- FIXME: all parameters are arbitrary at the moment
  $ and [ p > 0.01
        , p < 0.99
        , not $ isInfinite x0
        , not $ isInfinite x1
        ] ==> (if x0 < 1e6 then abs (x1 - x0) < 1e-6 else relativeError x1 x0 < 1e-12)
  where
    x0 = quantile      d (1 - p)
    x1 = complQuantile d p

-- Quantile is inverse of CDF
quantileIsInvCDF :: (Param d, ContDistr d) => T d -> d -> Double01 -> Property
quantileIsInvCDF _ d (Double01 p) =
  and [ p > m_tiny
      , p < 1
      , x > m_tiny
      , dens > 0
      ] ==>
    ( counterexample (printf "Quantile      = %g" x )
    $ counterexample (printf "Probability   = %g" p )
    $ counterexample (printf "Probability'  = %g" p')
    $ counterexample (printf "Rel. error    = %g" (relativeError p p'))
    $ counterexample (printf "Abs. error    = %e" (abs $ p - p'))
    $ counterexample (printf "Expected err. = %g" err)
    $ counterexample (printf "Distance      = %i" (ulpDistance p p'))
    $ counterexample (printf "Err/est       = %g" (fromIntegral (ulpDistance p p') / err))
    $ ulpDistance p p' <= round err
    )
  where
    -- Algorithm for error estimation is taken from here
    --
    -- http://sepulcarium.org/posts/2012-07-19-rounding_effect_on_inverse.html
    dens = density    d x
    err  = eps + eps' * abs (x / p) * dens
    --
    x    = quantile   d p
    p'   = cumulative d x
    (eps,eps') = prec_quantile_CDF d

-- Test that quantile fails if p<0 or p>1
quantileShouldFail :: (ContDistr d) => T d -> d -> Double -> Property
quantileShouldFail _ d p =
  p < 0 || p > 1 ==> QC.monadicIO $ do r <- QC.run $ E.catch
                                              (False <$ (return $! quantile d p))
                                              (\(_ :: E.SomeException) -> return True)
                                       QC.assert r

-- Probability is in [0,1] range
probSanityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Bool
probSanityCheck _ d x = p >= 0 && p <= 1
  where p = probability d x

-- Check that discrete CDF is correct
discreteCDFcorrect :: (DiscreteDistr d) => T d -> d -> Int -> Int -> Property
discreteCDFcorrect _ d a b
  = counterexample (printf "CDF   = %g" p1)
  $ counterexample (printf "Sum   = %g" p2)
  $ counterexample (printf "Delta = %g" (abs (p1 - p2)))
  $ abs (p1 - p2) < 3e-10
  -- Avoid too large differeneces. Otherwise there is to much to sum
  --
  -- Absolute difference is used guard againist precision loss when
  -- close values of CDF are subtracted
  where
    n  = min a b
    m  = n + (abs (a - b) `mod` 100)
    p1 = cumulative d (fromIntegral m + 0.5) - cumulative d (fromIntegral n - 0.5)
    p2 = sum $ map (probability d) [n .. m]

logProbabilityCheck :: (DiscreteDistr d) => T d -> d -> Int -> Property
logProbabilityCheck _ d x
  = counterexample (printf "probability    = %g" p)
  $ counterexample (printf "logProbability = %g" logP)
  $ counterexample (printf "log p          = %g" (log p))
  $ counterexample (printf "eps            = %g" (abs (logP - log p) / max (abs (log p)) (abs logP)))
  $ or [ p == 0     && logP == (-1/0)
       , p < 1e-308 && logP < 609
         -- To avoid problems with roundtripping error in case
         -- when density is computed as exponent of logDensity we
         -- accept either inequality
       ,  (ulpDistance (log p) logP <= 32)
       || (ulpDistance p (exp logP) <= 32)
       ]
  where
    p    = probability d x
    logP = logProbability d x

-- | Parameters for distribution testing. Some distribution require
--   relaxing parameters a bit
class Param a where
  -- | Whether quantileIsInvCDF is enabled
  quantileIsInvCDF_enabled :: T a -> Bool
  quantileIsInvCDF_enabled _ = True
  -- | Precision for 'quantileIsInvCDF' test
  prec_quantile_CDF :: a -> (Double,Double)
  prec_quantile_CDF _ = (16,16)
  -- |
  prec_discreteCDF :: a -> Double
  prec_discreteCDF _ = 32 * m_epsilon
  -- | Precision of CDF's complement
  prec_complementCDF :: a -> Double
  prec_complementCDF _ = 1e-14

instance Param StudentT where
  -- FIXME: disabled unless incompleteBeta troubles are sorted out
  quantileIsInvCDF_enabled _ = False
instance Param BetaDistribution where
  -- FIXME: See https://github.com/bos/statistics/issues/161 for details
  quantileIsInvCDF_enabled _ = False
instance Param FDistribution where
  -- FIXME: disabled unless incompleteBeta troubles are sorted out
  quantileIsInvCDF_enabled _ = False

instance Param ChiSquared where
  prec_quantile_CDF _ = (32,32)

instance Param BinomialDistribution where
  prec_discreteCDF _ = 1e-13
instance Param CauchyDistribution
instance Param DiscreteUniform
instance Param ExponentialDistribution
instance Param GammaDistribution
instance Param GeometricDistribution
instance Param GeometricDistribution0
instance Param HypergeometricDistribution
instance Param LaplaceDistribution
instance Param LognormalDistribution where
  prec_quantile_CDF _ = (64,64)
instance Param NormalDistribution
instance Param PoissonDistribution
instance Param UniformDistribution
instance Param WeibullDistribution
instance Param a => Param (LinearTransform a)

----------------------------------------------------------------
-- Unit tests
----------------------------------------------------------------

unitTests :: TestTree
unitTests = testGroup "Unit tests"
  [ testAssertion "density (gammaDistr 150 1/150) 1 == 4.883311" $
      4.883311418525483 =~ density (gammaDistr 150 (1/150)) 1
    -- Student-T
  , testStudentPDF 0.3  1.34  0.0648215  -- PDF
  , testStudentPDF 1    0.42  0.27058
  , testStudentPDF 4.4  0.33  0.352994
  , testStudentCDF 0.3  3.34  0.757146   -- CDF
  , testStudentCDF 1    0.42  0.626569
  , testStudentCDF 4.4  0.33  0.621739
    -- Student-T General
  , testStudentUnstandardizedPDF 0.3    1.2  4      0.45 0.0533456  -- PDF
  , testStudentUnstandardizedPDF 4.3  (-2.4) 3.22 (-0.6) 0.0971141
  , testStudentUnstandardizedPDF 3.8    0.22 7.62   0.14 0.0490523
  , testStudentUnstandardizedCDF 0.3    1.2  4      0.45 0.458035   -- CDF
  , testStudentUnstandardizedCDF 4.3  (-2.4) 3.22 (-0.6) 0.698001
  , testStudentUnstandardizedCDF 3.8    0.22 7.62   0.14 0.496076
    -- F-distribution
  , testFdistrPDF  1  3   3     (1/(6 * pi)) -- PDF
  , testFdistrPDF  2  2   1.2   0.206612
  , testFdistrPDF  10 12  8     0.000385613179281892790166
  , testFdistrCDF  1  3   3     0.81830988618379067153 -- CDF
  , testFdistrCDF  2  2   1.2   0.545455
  , testFdistrCDF  10 12  8     0.99935509863451408041
  ]
  where
    -- Student-T
    testStudentPDF ndf x exact
      = testAssertion (printf "density (studentT %f) %f ~ %f" ndf x exact)
      $ eq 1e-5  exact  (density (studentT ndf) x)
    testStudentCDF ndf x exact
      = testAssertion (printf "cumulative (studentT %f) %f ~ %f" ndf x exact)
      $ eq 1e-5  exact  (cumulative (studentT ndf) x)
    -- Student-T General
    testStudentUnstandardizedPDF ndf mu sigma x exact
      = testAssertion (printf "density (studentTUnstandardized %f %f %f) %f ~ %f" ndf mu sigma x exact)
      $ eq 1e-5  exact  (density (studentTUnstandardized ndf mu sigma) x)
    testStudentUnstandardizedCDF ndf mu sigma x exact
      = testAssertion (printf "cumulative (studentTUnstandardized %f %f %f) %f ~ %f" ndf mu sigma x exact)
      $ eq 1e-5  exact  (cumulative (studentTUnstandardized ndf mu sigma) x)
    -- F-distribution
    testFdistrPDF n m x exact
      = testAssertion (printf "density (fDistribution %i %i) %f ~ %f [got %f]" n m x exact d)
      $ eq 1e-5  exact d
      where d = density (fDistribution n m) x
    testFdistrCDF n m x exact
      = testAssertion (printf "cumulative (fDistribution %i %i) %f ~ %f [got %f]" n m x exact d)
      $ eq 1e-5  exact d
      where d = cumulative (fDistribution n m) x


text:
{-# LANGUAGE DeriveDataTypeable    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternSynonyms       #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeInType            #-}
{-# LANGUAGE UndecidableInstances  #-}
{-# LANGUAGE ViewPatterns          #-}

module Backprop.Learn.Model.Stochastic (
    DO(..)
  , StochFunc(..)
  , FixedStochFunc, pattern FSF, _fsfRunDeterm, _fsfRunStoch
  , rreLU
  , injectNoise, applyNoise
  , injectNoiseR, applyNoiseR
  ) where

import           Backprop.Learn.Model.Class
import           Backprop.Learn.Model.Function
import           Control.Monad.Primitive
import           Data.Bool
import           Data.Kind
import           Data.Typeable
import           GHC.TypeNats
import           Numeric.Backprop
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import qualified Data.Vector.Storable.Sized            as SVS
import qualified Statistics.Distribution               as Stat
import qualified System.Random.MWC                     as MWC
import qualified System.Random.MWC.Distributions       as MWC

-- | Dropout layer.  Parameterized by dropout percentage (should be between
-- 0 and 1).
--
-- 0 corresponds to no dropout, 1 corresponds to complete dropout of all
-- nodes every time.
newtype DO (n :: Nat) = DO { _doRate :: Double }
  deriving (Typeable)

instance KnownNat n => Learn (R n) (R n) (DO n) where
    runLearn (DO r) _ = stateless (constVar (realToFrac (1-r)) *)
    runLearnStoch (DO r) g _ = statelessM $ \x ->
        (x *) . constVar . vecR <$> SVS.replicateM (mask g)
      where
        mask = fmap (bool 1 0) . MWC.bernoulli r

-- | Represents a random-valued function, with a possible trainable
-- parameter.
--
-- Requires both a "deterministic" and a "stochastic" mode.  The
-- deterministic mode ideally should approximate some mean of the
-- stochastic mode.
data StochFunc :: Maybe Type -> Type -> Type -> Type where
    SF :: { _sfRunDeterm :: forall s. Reifies s W => Mayb (BVar s) p -> BVar s a -> BVar s b
          , _sfRunStoch
              :: forall m s. (PrimMonad m, Reifies s W)
              => MWC.Gen (PrimState m)
              -> Mayb (BVar s) p
              -> BVar s a
              -> m (BVar s b)
          }
       -> StochFunc p a b
  deriving (Typeable)

instance Learn a b (StochFunc p a b) where
    type LParamMaybe (StochFunc p a b) = p
    type LStateMaybe (StochFunc p a b) = 'Nothing

    runLearn SF{..} = stateless . _sfRunDeterm
    runLearnStoch SF{..} g = statelessM . _sfRunStoch g

-- | Convenient alias for a 'StochFunc' (random-valued function with both
-- deterministic and stochastic modes) with no trained parameters.
type FixedStochFunc = StochFunc 'Nothing

-- | Construct a 'FixedStochFunc'
pattern FSF :: (forall s. Reifies s W => BVar s a -> BVar s b)
            -> (forall m s. (PrimMonad m, Reifies s W) => MWC.Gen (PrimState m) -> BVar s a -> m (BVar s b))
            -> FixedStochFunc a b
pattern FSF { _fsfRunDeterm, _fsfRunStoch } <- (getFSF->(getWD->_fsfRunDeterm,getWS->_fsfRunStoch))
  where
    FSF d s = SF { _sfRunDeterm = const d
                 , _sfRunStoch  = const . s
                 }
{-# COMPLETE FSF #-}

newtype WrapDeterm a b = WD { getWD :: forall s. Reifies s W => BVar s a -> BVar s b }
newtype WrapStoch  a b = WS { getWS :: forall m s. (PrimMonad m, Reifies s W) => MWC.Gen (PrimState m) -> BVar s a -> m (BVar s b) }

getFSF :: FixedStochFunc a b -> (WrapDeterm a b, WrapStoch a b)
getFSF SF{..} = ( WD (_sfRunDeterm N_)
                , WS (`_sfRunStoch` N_)
                )

-- | Random leaky rectified linear unit
rreLU
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> FixedStochFunc (R n) (R n)
rreLU d = FSF { _fsfRunDeterm = vmap' (preLU v)
              , _fsfRunStoch  = \g x -> do
                  α <- vecR <$> SVS.replicateM (Stat.genContVar d g)
                  pure (zipWithVector preLU (constVar α) x)
              }
  where
    v :: BVar s Double
    v = constVar (Stat.mean d)

-- | Inject random noise.  Usually used between neural network layers, or
-- at the very beginning to pre-process input.
--
-- In non-stochastic mode, this adds the mean of the distribution.
injectNoise
    :: (Stat.ContGen d, Stat.Mean d, Fractional a)
    => d
    -> FixedStochFunc a a
injectNoise d = FSF { _fsfRunDeterm = (realToFrac (Stat.mean d) +)
                    , _fsfRunStoch  = \g x -> do
                        e <- Stat.genContVar d g
                        pure (realToFrac e + x)
                    }

-- | 'injectNoise' lifted to 'R'
injectNoiseR
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> FixedStochFunc (R n) (R n)
injectNoiseR d = FSF { _fsfRunDeterm = (realToFrac (Stat.mean d) +)
                     , _fsfRunStoch  = \g x -> do
                         e <- vecR <$> SVS.replicateM (Stat.genContVar d g)
                         pure (constVar e + x)
                     }

-- | Multply by random noise.  Can be used to implement dropout-like
-- behavior.
--
-- In non-stochastic mode, this scales by the mean of the distribution.
applyNoise
    :: (Stat.ContGen d, Stat.Mean d, Fractional a)
    => d
    -> FixedStochFunc a a
applyNoise d = FSF { _fsfRunDeterm = (realToFrac (Stat.mean d) *)
                   , _fsfRunStoch  = \g x -> do
                       e <- Stat.genContVar d g
                       pure (realToFrac e * x)
                   }

-- | 'applyNoise' lifted to 'R'
applyNoiseR
    :: (Stat.ContGen d, Stat.Mean d, KnownNat n)
    => d
    -> FixedStochFunc (R n) (R n)
applyNoiseR d = FSF { _fsfRunDeterm = (realToFrac (Stat.mean d) *)
                   , _fsfRunStoch  = \g x -> do
                       e <- vecR <$> SVS.replicateM (Stat.genContVar d g)
                       pure (constVar e * x)
                   }


text:
{-# LANGUAGE FlexibleContexts #-}
module NeuralNet (
    NeuralNet(),
    mkNeuralNet,
    trainNN,
    nnSize,
    nnCostFunction,
    feedForwardNN,
    randInitialNN,
    expandLogicalArray,
    flattenNN,
    reshapeNN,

    matrix,
    vector
  ) where

import           Control.Monad
import qualified Data.Vector as V
import           Numeric.GSL.Minimization
import           Numeric.LinearAlgebra
import           System.Random
import           Test.QuickCheck.Arbitrary (arbitrary)
import           Test.QuickCheck.Gen

import RandomMonad

data NeuralNet = NeuralNet
  { matrices :: V.Vector (Matrix Double)
  } deriving (Eq, Show)

mkNeuralNet :: [Matrix R] -> NeuralNet
mkNeuralNet = NeuralNet . V.fromList

trainNN :: Int -> R
        -> NeuralNet -> Matrix R -> Matrix R -> NeuralNet
trainNN maxIter lambda nn xs ys = reshapeNN dims solVec
  where
    dims = nnSize nn
    wrappedCF inVec = let (cost, grad) = nnCostFunction (reshapeNN dims inVec) xs ys lambda
                       in (cost, flattenNN grad)
    f = fst . wrappedCF
    df = snd . wrappedCF
    (solVec, optPath) = minimizeVD ConjugatePR 1E-3 maxIter 1E-2 1E-4 f df (flattenNN nn)

matDimsToNNDims :: [(Int,Int)] -> (Int,[Int],Int)
matDimsToNNDims ds = (inputCount, hiddenCounts, outputCount)
  where
    inputCount = snd (head ds) - 1
    hiddenCounts = fst <$> init ds
    outputCount = fst (last ds)

nnDimsToMatDims :: (Int,[Int],Int) -> [(Int,Int)]
nnDimsToMatDims (inputCount, hiddenCounts, outputCount)
  = go inputCount (hiddenCounts ++ [outputCount])
  where
    go _ [] = []
    go x (y:ys) = (y, x+1) : go y ys

nnSize :: NeuralNet -> (Int,[Int],Int)
nnSize nn = matDimsToNNDims (size <$> V.toList (matrices nn))

flattenNN :: NeuralNet -> Vector Double
flattenNN nn = vjoin (map flatten (V.toList (matrices nn)))

reshapeNN :: (Int,[Int],Int) -> Vector Double -> NeuralNet
reshapeNN dims xs = NeuralNet (V.fromList mats)
  where
    mats = zipWith reshape (map snd matDims) subVecs

    sv i c = subVector i c xs
    subVecs = zipWith sv idxs counts

    matDims = nnDimsToMatDims dims
    counts = map (\(r,c) -> r*c) matDims
    idxs = scanl (+) 0 counts

-- Returns the cost, and a neural net of gradients
nnCostFunction :: NeuralNet -> Matrix R -> Matrix R -> Double
               -> (Double, NeuralNet)
nnCostFunction nn xs ys lambda = (jCost, gradNN)
  where
    m = fromIntegral (rows xs)
    mats = matrices nn
    (activations, linearSums) = feedForward xs mats
    activationN = head activations

    jMatrix = (-ys * log activationN) - ((1-ys) * log (1-activationN))
    regVal = let drop1SqSum mat = sumElements ((dropColumns 1 mat) ** 2)
                 sqSum = V.sum (V.map drop1SqSum mats)
              in sqSum * lambda / 2
    jCost = (sumElements jMatrix + regVal) / m

    deltaN = activationN - ys
    -- We start with the Nth delta in our accumulator and work backwards
    -- We drop the first linearSum because that corresponds to layer N and we are
    -- working with the previous layer.
    deltas = reverse $ foldl computeDelta [deltaN] (tail (zip (V.toList mats) linearSums))
    computeDelta ds@(dNext:_) (theta, lSumPrev) = let dCurrent = dNext <> (dropColumns 1 theta) * sigGrad lSumPrev
                                                   in dCurrent:ds

    rawGrads = reverse $ zipWith (\del act -> (tr del <> (1 ||| act))) deltas (tail activations)
    -- we only regularize the non-bias terms
    regularizers = map (\mat -> scalar lambda * (0 ||| dropColumns 1 mat)) (V.toList mats)
    grads = zipWith (\g r -> (g + r) / scalar m) rawGrads regularizers
    gradNN = NeuralNet (V.fromList grads)

feedForwardNN :: NeuralNet -> Matrix R -> Matrix R
feedForwardNN nn xs = head (fst (feedForward xs (matrices nn)))

-- returns the final activation and progressive linear combinations (in reverse order)
feedForward :: Matrix R -> V.Vector (Matrix R) -> ([Matrix R], [Matrix R])
feedForward xs mats = V.foldl f ([xs], []) mats
  where
    f (as@(aPrev:_), ls) mat = let lCurrent = (1 ||| aPrev) <> tr mat
                                in ((sigmoid lCurrent):as, lCurrent:ls)

sigmoid x = 1 / (1 + exp(-x))
sigGrad x = (sigmoid x) * (1 - (sigmoid x))

randInitialNN :: (Int, [Int], Int) -> Rand NeuralNet
randInitialNN dims = do
  let matDims = nnDimsToMatDims dims
  mats <- mapM (uncurry randEpsilonMat) matDims
  return (NeuralNet (V.fromList mats))

randEpsilonMat :: Int -> Int -> Rand (Matrix R)
randEpsilonMat r c = do
  seed <- randRandom
  let mat = uniformSample seed r (replicate c (0,1))
      l_in = r - 1
      l_out = c
      epsilonInit = sqrt 6 / sqrt (fromIntegral (l_in * l_out))
  return ((mat * 2 * epsilonInit) - epsilonInit)

mkLogicalArray :: Num a => Int -> Int -> [a]
mkLogicalArray len pos = take len (replicate frontLen 0 ++ [1] ++ repeat 0)
  where
    frontLen = pos - 1

expandLogicalArray :: Int -> [Int] -> Matrix R
expandLogicalArray len is = matrix len (concatMap (mkLogicalArray len) is)

--------------------------------------------------------------------------------
-- QuickCheck Gen Instances
--------------------------------------------------------------------------------
positiveInt :: Gen Int
positiveInt = suchThat arbitrary (> 0)

nnDimsGen :: Gen (Int,[Int],Int)
nnDimsGen = (>**<) positiveInt (listOf positiveInt) positiveInt

nnGen :: Gen NeuralNet
nnGen = nnDimsGen >>= nnWithDimsGen

nnWithDimsGen :: (Int, [Int], Int) -> Gen NeuralNet
nnWithDimsGen dims = do
  let matDims = nnDimsToMatDims dims
  mats <- mapM (uncurry matrixGen) matDims
  return (NeuralNet (V.fromList mats))

matrixGen :: Int -> Int -> Gen (Matrix R)
matrixGen r c = do
  xs <- replicateM (r*c) arbitrary
  return (matrix c xs)

--------------------------------------------------------------------------------
-- Misc helpers
--------------------------------------------------------------------------------
-- Lifted from https://hackage.haskell.org/package/checkers-0.4.7/docs/src/Test-QuickCheck-Instances-Tuple.html#%3E%2A%3C
{- | Generates a 2-tuple using its arguments to generate the parts.
-}
(>*<) :: Gen a -> Gen b -> Gen (a,b)
x >*< y = liftM2 (,) x y

{- | Generates a 3-tuple using its arguments to generate the parts.
-}
(>**<) :: Gen a -> Gen b -> Gen c -> Gen (a,b,c)
(>**<) x y z = liftM3 (,,) x y z


text:
{-# LANGUAGE TypeOperators #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}

{-|
Module      : Numeric.LinearAlgebra.Combinators
Description : Compose matrices
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module Numeric.LinearAlgebra.Combinators
    ( blockMatrix
    )
where

import           GHC.TypeLits
import           Numeric.LinearAlgebra.Static

-- | Generate a block matrix from static matrices.
blockMatrix
    :: (KnownNat m, KnownNat n, KnownNat p, KnownNat q)
    => L m n             -- ^ Upper left matrix
    -> L m p             -- ^ Upper right matrix
    -> L q n             -- ^ Lower left matrix
    -> L q p             -- ^ Lower right matrix
    -> L (m + q) (n + p) -- ^ Resulting block matrix
blockMatrix a b c d = a ||| b === c ||| d