text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}

import           Control.Monad
import           Control.Monad.Random
import           Data.Constraint                 (Dict (..))
import           Data.List                       (foldl')
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.Prelude.List
import qualified Data.Vector.Storable            as V
import           GHC.TypeLits
import qualified Numeric.LinearAlgebra.Static    as SA
import           Options.Applicative
import           System.Random.MWC
import           System.Random.MWC.Distributions
import           Unsafe.Coerce                   (unsafeCoerce)

import           Grenade

-- | Choose a backend: 0 BLAS, 1 HMATRIX
#define HMATRIX 0

netSpec :: SpecNet
netSpec = specFullyConnected 2 400 |=> specRelu1D 400 |=> specFullyConnected 400 1 |=> specLogit1D 1 |=> specNil1D 1

-- netSpec :: SpecNet
-- netSpec = specFullyConnected 2 4000 |=> specRelu1D 4000 |=> specFullyConnected 4000 1 |=> specLogit1D 1 |=> specNil1D 1

-- netSpec :: SpecNet
-- netSpec =
--       specFullyConnected 2 4500
--   |=> specTanh1D 4500
--   |=> specFullyConnected 4500 750
--   |=> specRelu1D 750
--   |=> specFullyConnected 750 150
--   |=> specRelu1D 150
--   |=> specFullyConnected 150 150
--   |=> specRelu1D 150
--   |=> specFullyConnected 150 30
--   |=> specRelu1D 30
--   |=> specFullyConnected 30 20
--   |=> specRelu1D 20
--   |=> specFullyConnected 20 10
--   |=> specRelu1D 10
--   |=> specFullyConnected 10 1
--   |=> specLogit1D 1
--   |=> specNil1D 1

-- -- -- | The definition for a feed forward network using the dynamic module. Note the nested networks. This network clearly is over-engeneered for this example!
-- netSpec :: SpecNet
-- netSpec =
--  specFullyConnected 2 40 |=> specTanh1D 40 |=>
--  specDropout 400 0.95 Nothing |=>
--  netSpecInner |=>
--  specFullyConnected 20 30 |=> specRelu1D 30 |=>
--  specFullyConnected 30 20 |=> specRelu1D 20 |=>
--  specFullyConnected 20 10 |=> specRelu1D 10 |=>
--  specFullyConnected 10 1 |=> specLogit1D 1 |=> specNil1D 1
--  where netSpecInner = specFullyConnected 40 30 |=> specRelu1D 30 |=> specFullyConnected 30 20 |=> specReshape1D2D 20 (2, 10) |=> specReshape2D1D (2, 10) 20 |=> specNil1D 20

goals :: [((RealNum, RealNum), RealNum)]
goals = [((0.50, 0.50), 0.25), ((-0.50, -0.50), 0.25)]

-- | Specifications can be built using the following interface also. Here one does not have to carefully pay attention to match the dimension inputs and outputs as when using specification directly.
buildNetViaInterface :: IO SpecConcreteNetwork
buildNetViaInterface =
  buildModelWith (NetworkInitSettings UniformInit HMatrix (Just 10)) (DynamicBuildSetup { printResultingSpecification = False })  $
  inputLayer1D 2 >>                            -- 1. Every model has to start with an input dimension
  fullyConnected 10 >> dropout 0.89 >> relu >> -- 2. Layers are simply added as desired. The input of each layer is determined automatically,
  fullyConnected 20 >> relu >>                 --    whereas the output is specified as in `fullyConnected 20`. Thus, the layer empits 20 signals
  fullyConnected 4 >> tanhLayer >>             --    which arethe input of the next layer.
  networkLayer (
    inputLayer1D 4 >> fullyConnected 10 >> relu >> fullyConnected 4 >> sinusoid
    ) >>
  reshape (4, 1, 1) >>                         -- 3. Reshape is ignored if the input and output is the same
  fullyConnected 1                             -- 4. The output is simply the number of signals the alst layer emits.

netTrain :: forall shapes len1 len2 layers o .
     (SingI (Last shapes), KnownNat len1, KnownNat len2, Head shapes ~ 'D1 len1, Last shapes ~ 'D1 len2)
  => Network layers shapes
  -> Optimizer o
  -> Int
  -> IO (Network layers shapes)
netTrain net0 op n = do
#if HMATRIX
    inps <- replicateM n $ do
      s  <- getRandom
      return $ S1D (SA.randomVector s SA.Uniform * 2 - 1 :: SA.R len1)
#else
    inps <- replicateM n $ do
        s  <- getRandom
        return $ S1DV $ SA.extract (SA.randomVector s SA.Uniform * 2 - 1 :: SA.R len1)
#endif
    let mkVal :: S ('D1 len1) -> Rational
        mkVal (S1D v) = mkVal (S1DV $ SA.extract v)
        mkVal (S1DV v) = if any (v `inCircleV`) goals
                         then 1
                         else 0
#if HMATRIX
        constr = S1D . fromRational
#else
        constr = S1DV . fromRational
#endif
    let outs = map (constr . mkVal) inps
    let trained = foldl' trainEach net0 (zip inps outs)
    return trained

  where trainEach !network (i,o) = train op network i o

netScore :: (KnownNat len, Head shapes ~ 'D1 len, Last shapes ~ 'D1 1) => Network layers shapes -> IO ()
netScore network = do
    let testIns = [ [ (x,y)  | x <- [0..50] ]
                             | y <- [0..20] ]
#if HMATRIX
        outMat  = fmap (fmap (\(x,y) -> (render (x/25-1) (y/10-1) . normx) $ runNet network (S1D $ SA.vector [x / 25 - 1,y / 10 - 1]))) testIns
#else
        outMat  = fmap (fmap (\(x,y) -> (render (x/25-1) (y/10-1) . normx) $ runNet network (S1DV $ V.fromList [x / 25 - 1,y / 10 - 1]))) testIns
#endif
    putStrLn $ unlines outMat

  where
    render x y n'  | x == 0 && y == 0 = '+'
                   | y == 0 = '-'
                   | x == 0 = '|'
                   | n' <= 0.2  = ' '
                   | n' <= 0.4  = '.'
                   | n' <= 0.6  = '-'
                   | n' <= 0.8  = '='
                   | otherwise = '#'

normx :: S ('D1 1) -> RealNum
normx (S1D r)  = SA.mean r
normx (S1DV v) = V.sum v / fromIntegral (V.length v)

testValues :: forall len shapes layers . (KnownNat len, Head shapes ~ 'D1 len, Last shapes ~ 'D1 1) => Network layers shapes -> IO ()
testValues network = do
  let n = 1000
#if HMATRIX
  inps <- replicateM n $ do
      s  <- getRandom
      return $ S1D (SA.randomVector s SA.Uniform * 2 - 1 :: SA.R len)
#else
  inps <- replicateM n $ do
       s  <- getRandom
       return $ S1DV $ SA.extract (SA.randomVector s SA.Uniform * 2 - 1 :: SA.R len)
#endif

  let mkVal :: S ('D1 len1) -> Integer
      mkVal (S1D v) = mkVal (S1DV $ SA.extract v :: S ('D1 len))
      mkVal (S1DV v) = if any (v `inCircleV`) goals
                         then 1
                         else 0
  let outs :: [Integer]
      outs = map mkVal inps
  let ress = zip outs (map (round . normx . runNet network) inps)
      correct = length $ filter id $ map (uncurry (==)) ress
      incorrect = length $ filter id $ map (uncurry (/=)) ress
      falsePositives = length $ filter id $ map (uncurry (\shd nn -> shd == 0 && nn == 1)) ress
      falseNegatives = length $ filter id $ map (uncurry (\shd nn -> shd == 1 && nn == 0)) ress
  putStr $ show correct  ++ " | "
  putStr $ show incorrect ++ " | "
  putStr $ show falsePositives ++ " | "
  putStrLn $ show falseNegatives ++ " | "

inCircle :: KnownNat n => SA.R n -> (SA.R n, RealNum) -> Bool
v `inCircle` (o, r) = SA.norm_2 (v - o) <= r

inCircleV :: V.Vector RealNum -> ((RealNum, RealNum), RealNum) -> Bool
v `inCircleV` ((x,y), r) = norm (V.zipWith (-) v circle) <= r
  where
    circle = V.fromList [x,y]
    norm x = sqrt (V.sum $ V.map (^ (2 :: Int)) x)

data FeedForwardOpts = FeedForwardOpts Int (Optimizer 'Adam)

feedForward' :: Parser FeedForwardOpts
feedForward' =
  FeedForwardOpts <$> option auto (long "examples" <> short 'e' <> value 10000)
                  <*> (OptAdam
                       <$> option auto (long "alpha" <> short 'r' <> value 0.001)
                       <*> option auto (long "beta1" <> value 0.9)
                       <*> option auto (long "beta2" <> value 0.999)
                       <*> option auto (long "epsilon" <> value 1e-4)
                       <*> option auto (long "lambda" <> value 1e-3))

main :: IO ()
main = do

  FeedForwardOpts examples rate <- execParser (info (feedForward' <**> helper) idm)
  putStrLn "| Nr | Correct | Incorrect | FalsePositives | FalseNegatives |"
  putStrLn "--------------------------------------------------------------"
  let nr = 10 :: Int
  mapM_
    (\n -> do
       putStr $ "| " ++ show n ++ " | "
#if HMATRIX
       SpecConcreteNetwork1D1D (net0 :: Network layers shapes) <- networkFromSpecificationWith (NetworkInitSettings HeEtAl HMatrix Nothing) netSpec
#else
       SpecConcreteNetwork1D1D (net0 :: Network layers shapes) <- networkFromSpecificationWith (NetworkInitSettings HeEtAl BLAS Nothing) netSpec
#endif
       -- We need to specify the actual number of output nodes, as our functions requiere that!
       case (unsafeCoerce (Dict :: Dict ()) :: Dict (('D1 1) ~ Last shapes)) of
         Dict -> do
           net <- netTrain net0 rate examples
           unsafeCoerce $ -- only needed as GADTs are enabled, which disallowes the type to escape and thus prevents the type inference to work. The result is not needed anyways.
             testValues net)
    [1 .. nr]

  -- Features of dynamic networks:
  SpecConcreteNetwork1D1D (net' :: Network layers shapes) <- networkFromSpecificationWith (NetworkInitSettings HeEtAl BLAS (Just 50000)) netSpec
  net2 <- netTrain net' rate 0
  let spec' = networkToSpecification net2
  putStrLn "String represenation of the network specification: "
  print spec'
  let serializedSpec = encode spec'   -- only the specification (not the weights) are serialized here! The weights can be serialized using the networks serialize instance!
  let _ = encode net2                 -- E.g. like this.
  case decode serializedSpec of
    Left err -> print err
    Right spec2 -> do
      print spec2
      SpecConcreteNetwork1D1D (net3 :: Network layers3 shapes3) <- networkFromSpecificationWith (NetworkInitSettings HeEtAl BLAS (Just 50000)) spec2
      net4 <- foldM (\n _ -> netTrain n rate examples) net3 [(1 :: Int)..30]
      case (unsafeCoerce (Dict :: Dict ()) :: Dict (('D1 1) ~ Last shapes3)) of
        Dict -> netScore net4

  -- -- There is a also nice interface available also
  -- SpecConcreteNetwork1D1D newNet <- buildNetViaInterface
  -- print (networkToSpecification newNet)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
{-# LANGUAGE DeriveFunctor #-}
module Net where

import Prelude hiding (id,(.))

import Numeric.LinearAlgebra hiding (range)
import Data.Traversable
import System.Random
import Data.Array

import Debug.Trace
import Control.Category

--A function with a deriviative
--a and b are assumed to be numbers, and thus their own tangent/cotangent spaces
data DifferentiableFunction a b =
  DiffFunction { dfRun :: (a -> b)
               , dfDerivative :: (a -> (a -> b))
               }

type DF = DifferentiableFunction

instance Category DifferentiableFunction where
  id = DiffFunction id (const id)
  (DiffFunction f df) . (DiffFunction g dg) =
    DiffFunction (f . g) (\x -> (df (g x)) . (dg x))

simpleDiff :: (Num a) => DifferentiableFunction a b -> (a -> b)
simpleDiff df = (\x -> dfDerivative df x 1)

{-
instance Category DifferentiableFunction where
  id = DiffFunction id id
  (DiffFunction f df) . (DiffFunction g dg) = DiffFunction (f . g) ((f . dg)??df)
-}
data Layer c = Layer
  { layerWeights :: c
  , layerActivation :: DifferentiableFunction Double Double
  } deriving (Functor)

instance (Show c) => Show (Layer c) where
  show = show . layerWeights

data Net c = Net
  { netLayerSizes :: [Int]
  , netContent :: [Layer c]} deriving (Functor,Show)

netParameters :: Net c -> [c]
netParameters = (map layerWeights) . netContent

type NetShape = Net ()
type DoubleNet = Net (Matrix Double)

runLayer :: Layer (Matrix Double) -> Vector Double -> Vector Double
runLayer Layer{layerWeights=inputs,layerActivation=act} v =
  cmap (dfRun act) $ inputs #> v

runNet :: Net (Matrix Double) -> Vector Double -> Vector Double
runNet net input =
  --first layer needs a fake input that's always 1
  let realInput = vjoin [input,1.0]
  in foldl (flip runLayer) realInput (netContent net)

createSingle :: Int -> Int -> Double -> Vector Double
createSingle n i x = assoc n 0 [(i,x)]

runAndDiffLayer :: Layer (Matrix Double) -> Vector Double -> (Vector Double,Array (Int,Int) (Vector Double),Matrix Double)
runAndDiffLayer Layer{layerWeights=weights,layerActivation=act} input =
  let dact = simpleDiff act
      linearValues = weights #> input
      output = cmap (dfRun act) linearValues
      (nR,nC) = size weights
      limit = (nR-1,nC-1)

      deriveArrayValue r c = (dact $ atIndex linearValues r)*(atIndex input c)

      derivArray =
        listArray ((0,0),limit) $
        (\(r,c) -> createSingle nR r $ deriveArrayValue r c) <$> range ((0,0),limit)
      dactMatrix = (diag  (cmap dact linearValues))
  in (output,derivArray,dactMatrix <> weights)

runAndDiffNet :: Net (Matrix Double) -> Vector Double -> (Vector Double,[Array (Int,Int) (Vector Double)])
runAndDiffNet net input =
  let realInput = vjoin [input,1.0]
  in foldl (\(newIn,prevDiff) l ->
               let (nextIn,layerDiff,layerVDiff) = (runAndDiffLayer l newIn)
               in (nextIn,((fmap (layerVDiff #>)) <$> prevDiff) ++ [layerDiff]))
     (realInput,[]) $ netContent net

runAndDiffWithError :: Net (Matrix Double) -> Vector Double -> DifferentiableFunction (Vector Double) Double -> (Vector Double,Double,[Matrix Double])
runAndDiffWithError net input errFunc =
  let (out,grad) = runAndDiffNet net input
      derr = dfDerivative errFunc
  in (out,dfRun errFunc out,((fromArray2D . fmap (\v -> (derr out v))) <$> grad))

modifyNet :: Net a -> [b] -> Net b
modifyNet net contents =
  net{netContent = zipWith (\m l -> m <$ l) contents $ netContent net}

updateNet :: Net (Matrix Double) -> [Matrix Double] -> Double -> Net (Matrix Double)
updateNet net grad c =
  let param = netParameters net
      newParams = zipWith (\oldM gradM -> oldM + scale c gradM) param grad
  in modifyNet net newParams

tanhFunc :: DifferentiableFunction Double Double
tanhFunc = DiffFunction f (\y c -> c*(let z = f y in 1-z^2))
  where
    f :: Double -> Double
    f = tanh

--layerAsDiffFunc :: Layer (Matrix Double) -> DifferentiableFunction (Matrix Double,Vector Double) (Vector Double)

--netAsDiffFunction :: Net (Matrix Double) -> DifferentiableFunction ([Matrix Double],Vector Double) (Vector Double)

randomSimpleNet :: Int -> Int -> Int -> IO (Net (Matrix Double))
randomSimpleNet numIn numOut internal=
  do layer1 <- Layer <$> (scale (1/(fromIntegral (numIn+1))) <$> randn internal (numIn+1)) <*> pure tanhFunc

     layer2 <- Layer <$> (scale (1/(fromIntegral internal)) <$> randn numOut internal) <*> pure tanhFunc

     return $ Net [numIn,internal,numOut] [layer1,layer2]

ext: hs

max_stars_repo_name: jkelleyy/sinenet

max_issues_repo_name: jkelleyy/sinenet

max_forks_repo_name: jkelleyy/sinenet


text:
{-# LANGUAGE TemplateHaskell
           , TupleSections
           , DeriveGeneric
           , DeriveAnyClass
           , OverloadedLists #-}
module Physics.Collision where

import Numeric.LinearAlgebra (Vector, Matrix, (<>), scale, dot, cross, (#>))
import Control.Lens
import Numeric.LinearAlgebra.Data
import Data.Vector.Storable (init)
import GHC.Generics
import qualified Data.Vector as V
import qualified Data.Vector.Storable as VS
import Data.Function (on)

import Control.DeepSeq
import Utils
import Prelude hiding (init)

data Contact = Contact { _contactPoint :: Vector Float
                       , _contactNormal :: Vector Float
                       , _penetration :: Float
                       , _sign :: Float } deriving Show

makeLenses ''Contact

makeContactBasis cont =
  let ys = [0,1,0]
      y2 = [0,0,1]
      x  = cont^.contactNormal
      zs = normalize $ cross x ys
      z  = if VS.any isNaN zs
           then normalize $ cross x y2
           else zs
      y  = normalize $ cross z x
  in fromRows [x,y,z]

data CollisionData = CollisionData { _contacts :: [Contact]
                                   , _left :: Int }

makeLenses ''CollisionData

data Primitive = Sphere { _radius     :: Float
                        , _primoffset :: Matrix Float }
               | Plane  { _normal     :: Vector Float
                        , _planeoff   :: Float }
               | Box    { _points     :: V.Vector (Vector Float)
                        , _faces      :: V.Vector Squad }
                 deriving Show

data PrimitiveComputation = CSphere { _r    :: Float
                                    , _poff :: Vector Float }
                          | CPlane  { _cnormal   :: Vector Float
                                    , _cplaneoff :: Float}
                          | CBox    { _cpoints :: V.Vector (Vector Float)
                                    , _quads :: V.Vector Quad } deriving Show

data Quad = Quad { _qnorm   :: Vector Float
                 , _qoff    :: Float} deriving Show

quadToCPlane (Quad n o) = CPlane n o

data Squad = Squad { _snorm :: Vector Float
                   , _spoint:: Vector Float } deriving Show

sToQ t (Squad n p) = Quad newN o
  where
    newN = init . (t#>) . flip VS.snoc 0 $ n
    o    = dot newN . init . (t#>) . flip VS.snoc 1 $ p

preparePrimitive :: Matrix Float -> Primitive -> PrimitiveComputation
preparePrimitive t (Sphere r p1) = CSphere r ppos
  where ppos = init . (!!3) . toColumns $ t <> p1
preparePrimitive t (Box p q) = CBox (V.map adjust p) $ V.map (sToQ t) q
  where adjust = init . (t #>) . flip VS.snoc 1
preparePrimitive _ (Plane n o) = CPlane n o

tolerance = 0.05

quadDepth (Quad n o) p =
  (dot n p - o, n)

boxPtcoll (CBox _ qs) p cdata =
  let dots = V.map (flip quadDepth p) qs
      (d,n) = V.minimumBy (compare `on` fst) dots
  in if V.any ((>=0) . fst) dots
     then cdata
     else over contacts (Contact p n d 1:) cdata

quadPenetration (Quad n o) = planePenetration (CPlane n o)

planePenetration (CPlane n o) p cdata =
  let vertDist = dot p n
      cPoint   = p -- scale (vertDist - o) n + p
      contact  = Contact cPoint n (o - vertDist) 1
  in if cdata^.left > 0 && vertDist <= o + tolerance
     then over contacts (contact:) . over left (subtract 1) $ cdata
     else cdata

-- checkAxe :: Vector Float -> Vector Float -> V.Vector
checkAxe pos r ps ix =
  let vals = V.map (VS.!ix) ps
  in (pos VS.! ix + r) >= minimum vals &&
     (pos VS.! ix - r) <= maximum vals

checkConvexAxe ps1 ps2 ix =
  let vals1 = V.map (VS.! ix) ps1
      vals2 = V.map (VS.! ix) ps2
  in maximum vals1 >= minimum vals2 &&
     minimum vals1 <= maximum vals2

axeList :: [Int]
axeList = [0,1,2]

separatingAxes :: PrimitiveComputation -> PrimitiveComputation -> Bool
separatingAxes (CSphere r pos) (CBox ps _) = all (checkAxe pos r ps)     axeList
separatingAxes (CBox ps _) (CSphere r pos) = all (checkAxe pos r ps)     axeList
separatingAxes (CBox ps _) (CBox ps2 _)    = all (checkConvexAxe ps ps2) axeList

detectCollision prim1 prim2 =
  collisionDetector prim1 prim2 (CollisionData [] 8)

quadSphereCollision (CSphere r p) q@(Quad n o) =
  (dot n p - r - o, n)

quadContact (CSphere r p) (dist, n) =
  let contactPos = p - scale (dist + r) n
  in Contact contactPos n (-dist) 1

collisionDetector (CSphere r1 pos1) (CSphere r2 pos2) cdata =
  let mid     = pos1 - pos2
      midnorm = norm mid
  in if cdata^.left <= 0 || midnorm <= 0.0 || midnorm >= r1 + r2
     then cdata
     else let normal     = scale (1/midnorm) mid
              contactPos = pos1 + scale 0.5 mid
              contactPen = r1 + r2 - midnorm
              contact    = Contact contactPos normal contactPen 1
          in over contacts (contact:) . over left (subtract 1) $ cdata
collisionDetector p@(CPlane _ _) s@(CSphere _ _) cdata =
  set (contacts.traverse.sign) (-1) $
  collisionDetector s p cdata
collisionDetector (CSphere r pos1) (CPlane n p2) cdata =
  if cdata^.left <= 0 then cdata
  else let dist = dot n pos1 - r - p2
       in if dist >= 0
          then cdata
          else let contactPos = pos1 - scale (dist + r) n
                   contact    = Contact contactPos n (-dist) 1
               in over contacts (contact:) cdata
collisionDetector (CBox ps _) p@(CPlane _ _) cdata =
  V.foldr (planePenetration p) cdata ps
collisionDetector s@(CSphere _ _) b@(CBox _ _) cdata =
  set (contacts.traverse.sign) (-1) $
  collisionDetector b s cdata
collisionDetector (CBox _ q) s@(CSphere _ _) cdata =
  let cols = V.map (quadSphereCollision s) q
      best = V.maximumBy (compare `on` fst) cols
  in if V.any ((>=0) . fst) cols
     then cdata
     else over contacts (quadContact s best:) cdata
collisionDetector p@(CPlane _ _) b@(CBox _ _) cdata =
  set (contacts.traverse.sign) (-1) $
  collisionDetector b p cdata
collisionDetector b1@(CBox ps qs) b2@(CBox ps2 qs2) cdata =
  V.foldr (boxPtcoll b2) (V.foldr (boxPtcoll b1) cdata ps2) ps
collisionDetector _ _ cdata = cdata

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# LANGUAGE ScopedTypeVariables #-}

-- | Basic types used throughout the library.
--
--   Includes a probability type (based on sampling functions, with a tweak for
--   collected samples) and associated instances, as well as types used for
--   MCMC.

module Math.Probably.Types where

import Control.Applicative
import Control.Monad.State.Strict
import Numeric.LinearAlgebra
import System.Random.Mersenne.Pure64

type Seed = PureMT

data Prob a =
    Sampler { unSampler :: Seed -> (a, Seed) }
  | Samples [a]

instance Functor Prob where
  fmap f (Sampler sf) = Sampler $ \rs -> let (x,rs') = sf rs in (f x, rs')
  fmap f (Samples xs) = Samples $ map f xs

instance Applicative Prob where
  pure x = Sampler (\rs-> (x, rs))
  (Sampler sff) <*> (Sampler sfx) = Sampler $ \rs ->
    let (f ,rs') = sff rs
        (x, rs'') = sfx rs'
    in (f x, rs'')
  _ <*> _ = error "Prob (<*>): unsupported pattern"

instance Monad Prob where
  return = pure
  (Sampler sf) >>= f = Sampler $ \rs->
    let (x, rs'::Seed) = sf rs
        nextProb = f x
    in case nextProb of
         Sampler g -> g rs'
         Samples xs -> primOneOf xs rs'

  (Samples xs) >>= f = Sampler $ \rs->
    let (x, rs'::Seed) = primOneOf xs rs
        nextProb = f x
    in case nextProb of
         Sampler g -> g rs'
         Samples ys -> primOneOf ys rs'

primOneOf :: [a] -> Seed -> (a, Seed)
primOneOf xs seed
  = let (u, nextSeed) = randomDouble seed
        idx = floor $ realToFrac u * realToFrac (length xs )
    in (xs !! idx, nextSeed)

type DiscreteParams   = Vector Int
type ContinuousParams = Vector Double
type Parameters       = (DiscreteParams, ContinuousParams)

type LogObjective     = Parameters -> Double
type Gradient         = ContinuousParams -> ContinuousParams

type Particle         = (ContinuousParams, ContinuousParams)

type Transition t     = StateT (Chain t) Prob Parameters

-- | State of a Markov chain.  Note that the objective function itself is
--   included in the state, which allows the possibility of annealing schedules
--   and the like.
data Chain t = Chain {
    parameterSpacePosition :: Parameters
  , objectiveFunction      :: Target
  , objectiveValue         :: Double
  , tunables               :: t
  }

-- | A target to sample, consisting of a log objective function and possibly
--   its gradient.  The gradient is only taken with respect to continuous
--   parameters.
data Target = Target {
    logObjective :: Parameters -> Double
  , gradient     :: Maybe (ContinuousParams -> ContinuousParams)
  }

-- | Convenience constructor for targets and gradients.
createTargetWithGradient :: LogObjective -> Gradient -> Target
createTargetWithGradient f g = Target f (Just g)

-- | Convenience constructor for targets without gradients.
createTargetWithoutGradient :: LogObjective -> Target
createTargetWithoutGradient f = Target f Nothing

handleGradient :: Maybe t -> t
handleGradient Nothing  = error "handleGradient: no gradient provided"
handleGradient (Just g) = g

-- | The dual-averaging implementation of NUTS requires quite a few tuning
--   parameters, which are all held in this type.
data DualAveragingParameters = DualAveragingParameters {
    mAdapt    :: !Int
  , delta     :: !Double
  , mu        :: !Double
  , gammaP    :: !Double
  , tau0      :: !Double
  , kappa     :: !Double
  , daStep    :: !Double
  , daStepAvg :: !Double
  , daH       :: !Double
  } deriving (Eq, Show)

defaultDualAveragingParameters :: Double -> Int -> DualAveragingParameters
defaultDualAveragingParameters e burnInPeriod = DualAveragingParameters {
    mu        = log (10 * e)
  , delta     = 0.5
  , mAdapt    = burnInPeriod
  , gammaP    = 0.05
  , tau0      = 10
  , kappa     = 0.75
  , daStep    = e
  , daStepAvg = e
  , daH       = 0
  }

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
module Main where
import           Control.Monad.State
import           Data.List           (sortBy, (\\))
import           Data.Matrix         ((<->), (<|>))
import qualified Data.Matrix         as M
import           Data.Monoid         ((<>))
import           Data.Ord            (comparing)
import           Data.Random         (sampleState, shuffleNofM)
import           Data.Text           (Text, unpack)
import qualified Data.Text.IO        as T (readFile)
import qualified Data.Vector         as V
import           Statistics.Sample   (mean, stdDev)
import           System.Environment  (getArgs)
import           System.IO           (hFlush, stdout)
import           System.Random       (mkStdGen)

{- massage the big string into a good dataset
-}

parse :: String -> M.Matrix Double
parse =
    M.fromLists . map (map read . words) . lines

{- Split the initial matrix into a vector of labels and a matrix of features
-}

splitData :: M.Matrix Double -> (V.Vector Double, M.Matrix Double)
splitData m =
    ( M.getCol 1 m
    , M.submatrix 1 (M.nrows m) 2 (M.ncols m) m
    )

prompt :: IO (String, String) -- The filepath and the algorithm number
prompt = do
    args <- getArgs
    putStrLn "Welcome to Nick Lawler's Feature Selection Algorithm"
    case args of
        [filename, algorithm] -> do
            putStrLn "Filename and Algorithm provided at command line\n"
            return (filename,algorithm)
        _ -> do
            putStr "Type in the name of the file to test : "
            hFlush stdout
            filename <- getLine
            putStrLn ""
            putStrLn "Type the number of the algorithm you want to run."
            putStrLn ""
            putStrLn "        1) Forward Selection"
            putStrLn "        2) Backward Elimination"
            putStrLn "        3) Nick’s Special Algorithm."
            putStr "                                    "
            hFlush stdout
            algorithm <- getLine
            return (filename,algorithm)

{- Z normalize the matrix one column at a time
-}

normalize :: M.Matrix Double -> M.Matrix Double
normalize =
    let normalizeCol col = V.map (\x -> (x - mean col) / stdDev col ) col
    in V.foldl1' (<|>) . V.map (M.colVector . normalizeCol) . toCols

main :: IO ()
main = do
    (filename, algorithm) <- prompt
    raw <- T.readFile filename
    let (labels, rawDataset) = splitData (parse (unpack raw))
    putStrLn $
          "This dataset has "
          ++ show (M.ncols rawDataset)
          ++ " features (not including the class attribute) with "
          ++ show (M.nrows rawDataset)
          ++ " instances.\n"
    putStr "Please wait while I normalize the data... "
    hFlush stdout
    let dataset = M.forceMatrix (normalize rawDataset)
    putStrLn "Done!"
    let algo = case algorithm of
                    "1" -> forwardSelection
                    "2" -> backwardElimination
                    "3" -> specialAlgorithm
                    _ -> error "Choose a proper algorithm!"
    putStrLn $ "Running nearest neighbor with all "
            ++ show (M.ncols dataset)
            ++ " features, using “leaving-one-out” evaluation, I get an accuracy of "
            ++ show (leaveOneOutCV [1.. M.ncols dataset] labels dataset)
            ++ "%\n"
    putStrLn "Beginning Search\n"
    (accuracy,featureSet) <- algo labels dataset
    putStrLn $ "Finished Search!! The best feature subset is "
            ++ show featureSet
            ++ ", which has an accuracy of "
            ++ show accuracy
            ++ "%"

{- The full cross validation pipeline

   After computing the nearest neighbor for each validation run, pair them with
   the true labels, count the number of times you computed the correct one and
   return the proportion of correct classifications
-}

leaveOneOutCV :: [Int]           -- features to consider
              -> V.Vector Double -- labels
              -> M.Matrix Double -- dataset
              -> Double          -- the percentage of computed labels that matched the true lable
leaveOneOutCV features labels dataset =
    let correctlyLabled = V.length
                        $ V.filter (== True)
                        $ V.zipWith (==) labels
                        $ allNearestNeighbors labels (getCols (V.fromList features) dataset)
    in 100 * (fromIntegral correctlyLabled / fromIntegral (M.nrows dataset))

{- Picks out the columns of 'matrix' specified by their 1 based indices in the array cols
-}

getCols :: V.Vector Int -> M.Matrix a -> M.Matrix a
getCols cols matrix =
    if   V.null cols then M.fromLists [[]]
    else V.foldl1' (<|>) $ V.map (M.colVector . (`M.getCol` matrix)) cols

{- Apply nearest neighbor to all (true label, example, training set) permutations
-}

allNearestNeighbors :: V.Vector Double -> M.Matrix Double -> V.Vector Double
allNearestNeighbors labels dataset =
    V.zipWith3 nearestNeighbor (leaveOneOutVec labels) (toRows dataset) (leaveOneOut dataset)

{- Returns a vector containing all possible removals of a single element from
   the argument vector
-}

leaveOneOutVec :: V.Vector a -> V.Vector (V.Vector a)
leaveOneOutVec vec =
    V.map (removeElement vec) (V.enumFromN 0 (V.length vec))

{- Returns a vector equivalent to the argument matrix except that 0-based
   element i has been removed
-}

removeElement :: V.Vector a -> Int -> V.Vector a
removeElement vec i =
    let (top, bottom) = V.splitAt i vec
    in  top <> V.slice 1 (V.length bottom - 1) bottom

{- Returns a vector containing all possible removals of a single row from the
   argument matrix
-}

leaveOneOut :: M.Matrix a -> V.Vector (M.Matrix a)
leaveOneOut matrix =
    V.map (removeRow matrix) (V.enumFromN 1 (M.nrows matrix))

{- Returns a copy of the argument matrix missing the 1-indexed row provided as
   the second argumetn
-}

removeRow :: M.Matrix a -> Int -> M.Matrix a
removeRow matrix rowNum
    | rowNum == 1              = rowsBelow
    | rowNum == M.nrows matrix = rowsAbove
    | otherwise                = rowsAbove <-> rowsBelow
    where rowsAbove = M.submatrix 1 (rowNum - 1) 1 (M.ncols matrix) matrix
          rowsBelow = M.submatrix (rowNum+1) (M.nrows matrix) 1 (M.ncols matrix) matrix

{- The core nearest neighbor algorithm.

   After turning the dataset into a vector of rows, it calculates the euclidian
   distance between each row and the example, pairs that distance with the
   provided label, then extracts the label of the minimum-scoring example

-}
nearestNeighbor :: V.Vector Double -- the labels
                -> V.Vector Double -- the example to be classified
                -> M.Matrix Double -- the dataset
                -> Double          -- the computed label
nearestNeighbor labels example =
      fst
    . V.minimumBy (comparing snd)
    . V.zip labels
    . V.map (euclidianDistance example)
    . toRows

{- Calculate the euclidian distance between two numeric vectors
-}

euclidianDistance :: Num a => V.Vector a -> V.Vector a -> a
euclidianDistance v1 v2 =
      V.sum (V.map (^(2 :: Int)) (V.zipWith (-) v1 v2))

{- Represent a matrix as a vector of row vectors
-}

toRows :: M.Matrix a -> V.Vector (V.Vector a)
toRows matrix =
    V.map (`M.getRow` matrix) (V.enumFromN 1 (M.nrows matrix))

{- Represent a matrix as a vector of Column vectors
-}

toCols :: M.Matrix a -> V.Vector (V.Vector a)
toCols = toRows . M.transpose

{- It turns out that the vast majority of the code for the forward and backward
   searches is essentially the same. We can reuse the main "featureSearch" section
   by abstracting out a function to choose which features are examined in the
   inner loop, another function to determine whether we add or remove a feature
   on each iteration, and the inital set of features.
-}

featureSearch :: ([Int] -> [Int] -> [Int]) -- how to choose which features to consider
              -> ([Int] -> [Int] -> [Int]) -- how to add the next feature
              -> [Int]                     -- starting features
              -> V.Vector Double
              -> M.Matrix Double
              -> IO (Double, [Int])
featureSearch choose add initial labels dataset = do
    (_,_,acc,features) <- execStateT computeFeatures (50,initial,50,initial)
    return (acc,features)
    where
        -- Haskell forces us to be clear about our intention to use 4 mutable variables
        -- even then, its only a functional simulation of mutable state
        computeFeatures :: StateT (Double, [Int], Double, [Int]) IO ()
        computeFeatures =
            forM_ [1.. M.ncols dataset] $ \_ -> do
                (lastBestAcc, lastBestFeatures, bestAcc, bestFeatures) <- get
                let featuresToConsider = [1.. M.ncols dataset] `choose` lastBestFeatures
                accSets <- forM featuresToConsider $ \k -> do
                            let candidateSet = lastBestFeatures `add` [k]
                                (accuracy,set) = case candidateSet of
                                                [] -> (leaveOneOutCV [1..M.ncols dataset] labels dataset,[1..M.ncols dataset])
                                                _  -> (leaveOneOutCV candidateSet labels dataset,candidateSet)
                            putStrLnM ("       Using feature(s) " ++ show set ++ " accuracy is " ++ show accuracy ++ "%")
                            return (accuracy, set)
                let (newAcc, newFeatures) = maximum accSets -- maximum looks at left entry of tuple
                putStrLnM  ""
                when (newAcc < lastBestAcc) $ -- possible corner case when equal
                    putStrLnM "(Warning, Accuracy has decreased! Continuing search in case of local maxima)"
                putStrLnM ("Feature set " ++ show newFeatures ++ " was best, accuracy " ++ show newAcc ++ "%\n")
                put (newAcc, newFeatures, max newAcc bestAcc, if newAcc > bestAcc then newFeatures else bestFeatures)

-- to choose the features considered, we remove the ones we've already looked at
forwardSelection :: V.Vector Double -> M.Matrix Double -> IO (Double, [Int])
forwardSelection = featureSearch (\\) (++) []

-- we're only interested in the features that are left after the last iteration
-- so we use `seq`, which always returns its second argument.
backwardElimination :: V.Vector Double -> M.Matrix Double -> IO (Double,[Int])
backwardElimination labels dataset = featureSearch seq (\\) [1..M.ncols dataset] labels dataset

{- in each inner loop, try out the top 3 scoring feature subsets on the full dataset.
   pick the one with the highest score and maintain that score going forward
-}
specialAlgorithm :: V.Vector Double -> M.Matrix Double -> IO (Double,[Int])
specialAlgorithm labels dataset = do
    putStrLn "Taking a random sample of the original dataset"
    let sampleSize = M.nrows dataset `div` 10
    (_,_,acc,features) <- execStateT  (uncurry computeFeatures (randomSample sampleSize labels dataset)) (50,[],50,[])
    return (acc,features)
    where
        computeFeatures :: V.Vector Double -> M.Matrix Double -> StateT (Double, [Int], Double, [Int]) IO ()
        computeFeatures labels' dataset' =
            forM_ [1.. M.ncols dataset'] $ \_ -> do
                (lastBestAcc, lastBestFeatures, bestAcc, bestFeatures) <- get
                let featuresToConsider = [1.. M.ncols dataset'] \\ lastBestFeatures
                accSets <- forM featuresToConsider $ \k -> do
                            let candidateSet = lastBestFeatures ++ [k]
                                (accuracy,set) = case candidateSet of
                                                [] -> (leaveOneOutCV [1..M.ncols dataset'] labels' dataset',[1..M.ncols dataset'])
                                                _  -> (leaveOneOutCV candidateSet labels' dataset',candidateSet)
                            putStrLnM ("       Using feature(s) " ++ show set ++ " accuracy is " ++ show accuracy ++ "%")
                            return (accuracy, set)
                let accSets' = topNOnFullDataset 3 labels dataset accSets
                    (newAcc, newFeatures) = maximum accSets' -- maximum looks at left entry of tuple
                putStrLnM  ""
                when (newAcc < lastBestAcc) $ -- possible corner case when equal
                    putStrLnM "(Warning, Accuracy has decreased! Continuing search in case of local maxima)"
                putStrLnM ("Feature set " ++ show newFeatures ++ " was best, full dataset accuracy " ++ show newAcc ++ "%\n")
                put (newAcc, newFeatures, max newAcc bestAcc, if newAcc > bestAcc then newFeatures else bestFeatures)

{- Compute the full-dataset validation score of the features subsets that scored
   the best on the reduced-size dataset
-}

topNOnFullDataset :: Int -> V.Vector Double -> M.Matrix Double -> [(Double,[Int])] -> [(Double,[Int])]
topNOnFullDataset n labels dataset =
    let scoreWithFeatures labels' dataset' (_,features) = (leaveOneOutCV features labels' dataset', features)
    in map (scoreWithFeatures labels dataset) . take n . sortBy (flip compare)

{- Return a random n-length subset of the supplied label vector and dataset matrix
-}

randomSample :: Int -> V.Vector Double -> M.Matrix Double -> (V.Vector Double, M.Matrix Double)
randomSample n labels dataset =
    let (labels',dataset') = unzip . fst $ sampleState (shuffleNofM n (M.nrows dataset) (zip (V.toList labels) (M.toLists dataset))) (mkStdGen 366)
    in (V.fromList labels', M.fromLists dataset')

-- quick dataset access in repl
small34 :: IO Text
small34 = T.readFile "data/cs_170_small34.txt"

large34 :: IO Text
large34 = T.readFile "data/cs_170_large34.txt"

small80 :: IO Text
small80 = T.readFile "data/cs_170_small80.txt"

large80 :: IO Text
large80 = T.readFile "data/cs_170_large34.txt"

-- easy printing in IO Monads
putStrLnM :: MonadIO m => String -> m ()
putStrLnM = liftIO . putStrLn

max_stars_repo_name: nicklawls/features

max_issues_repo_name: nicklawls/features

max_forks_repo_name: nicklawls/features


text:
{-# LANGUAGE LambdaCase #-}
module Statistics.Graph
    ( graphF
    , writeGraph
    ) where

import           Universum hiding (unlines)

import           Control.Foldl (Fold (..))
import           Data.Graph.Inductive.Graph (Graph (mkGraph))
import           Data.Graph.Inductive.PatriciaTree (Gr)
import           Data.GraphViz (DotGraph)
import qualified Data.GraphViz as G
import qualified Data.GraphViz.Attributes.Complete as A
import           Data.GraphViz.Commands.IO (hPutDot)
import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as M
import           Data.Set (Set)
import qualified Data.Set as S
import           Pos.Util (withTempFile)
import           Prelude (unlines)
import           System.Exit (ExitCode (ExitSuccess))
import           System.IO (hGetContents, hPutStrLn)
import           System.Process (readProcessWithExitCode)

import           JSONLog (IndexedJLTimedEvent)
import           Statistics.Block (BlockHeader (..), blockChain, blockHeadersF)
import           Types

graphF :: Fold IndexedJLTimedEvent (DotGraph Int)
graphF = f <$> blockHeadersF
  where
    f :: Map BlockHash BlockHeader -> DotGraph Int
    f m =
        let nodes = zip [1..] $ [bh | (_, bh) <- M.toList m]
            h2i   = M.fromList [(bhHash bh, i) | (i, bh) <- nodes]
            edges = do
                (i, bh) <- nodes
                let h = bhHash bh
                    h' = bhPrevBlock $ m M.! h
                guard $ M.member h' h2i
                return (i, h2i M.! h', ())
            g     = mkGraph nodes edges :: Gr BlockHeader ()
        in  G.graphToDot G.nonClusteredParams { G.fmtNode = fmt } g
      where
        chain :: Set BlockHash
        chain = blockChain m

        fmt :: (Int, BlockHeader) -> G.Attributes
        fmt (_, bh) = [ A.Label $ A.StrLabel $ toLText $ unlines labels
                      , A.FillColor (if S.member (bhHash bh) chain
                            then [A.WC (A.X11Color G.Yellow) Nothing]
                            else [A.WC (A.X11Color G.LightGray) Nothing])
                      , A.Style [A.SItem A.Filled []]
                      ]
          where
            labels :: [String]
            labels = [ '#' : take 6 (toString $ bhHash bh)
                     , show (bhNode bh)
                     , show (bhSlot bh)
                     ]

writeGraph :: FilePath -> DotGraph Int -> IO Bool
writeGraph f g = withTempFile "." "graph.dot" $ \_ h -> do
    hPutDot h g
    b <- G.isGraphvizInstalled
    case b of
        True -> do
            input <- hGetContents h
            ex <- view _1 <$> readProcessWithExitCode "dot" ["-Tpng", "-o" ++ f] input
            case ex of
                ExitSuccess -> return True
                _           -> hPutStrLn Universum.stderr ("Creating the graph failed, " ++ show ex) >> return False
        False -> hPutStrLn stderr "Cannot produce graph without dot. Please install graphviz." >> return False

max_stars_repo_name: Davdegreat/cardano-sl

max_issues_repo_name: Davdegreat/cardano-sl

max_forks_repo_name: Davdegreat/cardano-sl


text:
-- Copyright 2016 TensorFlow authors.
--
-- Licensed under the Apache License, Version 2.0 (the "License");
-- you may not use this file except in compliance with the License.
-- You may obtain a copy of the License at
--
--     http://www.apache.org/licenses/LICENSE-2.0
--
-- Unless required by applicable law or agreed to in writing, software
-- distributed under the License is distributed on an "AS IS" BASIS,
-- WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-- See the License for the specific language governing permissions and
-- limitations under the License.

{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeFamilies #-}
{- | Rendering of TensorFlow operations as Haskell functions.

The basic type signature generated for each op is:

> {constraints} => {mandatory attrs} -> {input tensors} -> {output tensors}

where:

* @{mandatory attrs}@ is of the form @A_1 -> ... -> A_N@, where each @A@ is an
 op attribute that doesn't have a default and can't be inferred from other
 inputs.

* @{constraints}@ restrict the type parameters of the input and output tensors
 (for example: 'TensorType' or 'OneOf').

* @{input tensors}@ is of the form @T_1 -> ... -> T_N@, where each @T@ is of
the form @Tensor Ref a@ or @Tensor v a@ (or a list of one of those types),
and @a@ is either a concrete type or a (constrained) type variable.

* @{output tensors}@ is of the form @(T_1,...,T_N)@ for "pure" ops, and
@Build (T_1,...,T_N)@ for "stateful" ops.  An op is considered "stateful" if
it takes a @Tensor Ref@ or @Tensor v ResourceHandle@ as input, or if it's
explicitly marked \"Stateful\" in its @REGISTER_OP@ definition.  (If there
are no outputs, it is either @ControlNode@ or @Build ControlNode@.)
-}

module TensorFlow.OpGen
  ( OpGenFlags(..)
  , docOpList
  , flagParser)
  where

import Data.Foldable (toList)
import Data.Maybe (fromMaybe)
import Data.ProtoLens (def, showMessage)
import Data.List (sortOn)
import Data.List.NonEmpty (NonEmpty)
import qualified Data.List.NonEmpty as NE
import Lens.Family2 ((^.), (.~), (&), view)
import Options.Applicative (Parser, help, long, strOption, value)
import Proto.Tensorflow.Core.Framework.OpDef
  ( OpList
  , OpDef
  , attr
  , inputArg
  , name
  , op
  , outputArg
  )
import Proto.Tensorflow.Core.Framework.Types (DataType(..))
import System.FilePath (takeBaseName)
import TensorFlow.OpGen.ParsedOp
import Text.PrettyPrint.Mainland
  ( Doc
  , (<>)
  , (<+>)
  , (</>)
  , (<+/>)
  , brackets
  , comma
  , commasep
  , dquotes
  , empty
  , enclose
  , flatten
  , folddoc
  , hang
  , indent
  , parens
  , sep
  , stack
  , strictText
  , tuple
  )
import qualified Data.Set as Set
import qualified Data.Text as Text

data OpGenFlags = OpGenFlags
     { outputFile :: String
     , prefix :: String
     , excludeList :: String
     }

flagParser :: Parser OpGenFlags
flagParser = OpGenFlags
     <$> strOption (mconcat [ long "output"
                            , help "File to write."
                            ])
     <*> strOption (mconcat [ long "prefix"
                            , help "Haskell package prefix to use"
                            ])
     <*> strOption (mconcat [ long "exclude_list"
                            , value ""
                            , help "Comma separated Ops names to ignore"
                            ])

docOpList :: OpGenFlags -> OpList -> Doc
docOpList flags opList =
  stack [ "{-# LANGUAGE ConstraintKinds #-}"
        , "{-# LANGUAGE DataKinds #-}"
        , "{-# LANGUAGE FlexibleContexts #-}"
        , "{-# LANGUAGE FlexibleInstances #-}"
        , "{-# LANGUAGE OverloadedStrings #-}"
        , "{-# LANGUAGE ScopedTypeVariables #-}"
          -- Avoids reports about shadowing standard library names.
        , "{-# OPTIONS_GHC -fno-warn-name-shadowing #-}"
          -- eqLengthGuard never returns false and dies instead.
        , "{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}"
        , "module" <+> strictText moduleName <+> "where"
        , empty
        , imports
        , empty
        , folddoc (\x y -> x </> empty </> y)
                  (map renderOpAndExtras $
                   sortOn (view name) $
                   filter (not . flip elem exclusions . view name) $
                   toList $ opList ^. op)
        ]
  where moduleName =
            Text.pack (prefix flags) <> "." <> camelCase
             -- Discards the optional trailing _ops_op_lib
            (fromMaybe shortName (Text.stripSuffix "_ops_op_lib" shortName))
        shortName = Text.pack (takeBaseName $ outputFile flags)
        exclusions = Text.splitOn "," $ Text.pack $ excludeList flags
        renderOpAndExtras o = renderOp (parseOp o) </> extras o

imports :: Doc
imports = stack [
      "import Data.ByteString (ByteString)"
    , "import Data.Complex (Complex)"
    , "import Data.Int (Int8, Int16, Int32, Int64)"
    , "import Data.Proxy (Proxy(Proxy))"
    , "import Data.Word (Word8, Word16, Word32, Word64)"
    , "import Lens.Family2 ((.~), (&))"
    , "import TensorFlow.Build"
    , "import TensorFlow.BuildOp"
    , "import TensorFlow.Tensor"
    , "import TensorFlow.Types"
    ]

renderHaskellName, renderTFName, renderQuotedTFName :: Name -> Doc
renderHaskellName = strictText . unHaskellName . haskellName
renderTFName = strictText . unTFName . tfName
renderQuotedTFName = dquotes . renderTFName

-- | Generate the source code for a single op.
-- For example:
--
-- -- | {haddock comment}
-- foo :: {type sig}
-- foo attr1 attr2 input1 input2 | eqLengthGuard [...] = {function body}
renderOp :: ParsedOp -> Doc
renderOp pOp = stack $
    [ haddocks
    -- Prevent unreasonably long compilation times on ghc-7.10, due
    -- to stack calling "-dump-hi" which (unnecessarily) includes the
    -- inlining information, and is large for ops with many arguments.
#if __GLASGOW_HASKELL__ < 800
    , "{-# NOINLINE" <+> n <+> "#-}"
#endif
    , n <+> "::" <+> hang 0 (typeSig empty pOp)
    , n <+> "=" <+> n <> "' id"
    , n' <+> "::" <+> hang 0 (typeSig "OpParams ->" pOp)
    , n' <+> hang 0 args <+> "|" <+> funcGuard listSizeAttrs
                <+> "=" </>  -- args are indented
                    -- the body needs to be indented wrt the name
                    indent indentation (functionBody pOp)
    ] ++ whereClause listSizeAttrs
  where
    n = renderHaskellName $ parsedOpName pOp
    n' = n <> "'"
    listSizeAttrs = inferredListSizeAttrs pOp
    args = sep $ "op'options"
               : (map renderHaskellName
                    $ map attrName (explicitInputAttrs pOp)
                    ++ map parsedArgName (parsedInputs pOp))
    haddocks = "-- |" <+> multilineComment (parsedOpSummary pOp) (parsedOpDescription pOp)

-- | A check that all lists of the given size have the given length.
-- For example:
--   eqLengthGuard [("N", [("input1", length input1), ("input2", length input2)])]
funcGuard :: [Attr (NonEmpty Name)] -> Doc
funcGuard attrs = "eqLengthGuard" <+> brackets (commasep entries)
      where
        entries =
            [ parens $ nAttr <> comma <+>
              brackets (commasep $ toList $
                            map renderTensorName (toList $ attrInfo a))
            | a <- attrs
            , let nAttr = renderQuotedTFName (attrName a)
            ]
        renderTensorName x = parens $ renderQuotedTFName x <> comma <+>
                        "length" <+> renderHaskellName x

-- | Define the implicit list length attributes.
-- For example:
--   where
--     n1 = fromIntegral (length input1) :: Int64
--     n2 = fromIntegral (length input2) :: Int64
whereClause :: [Attr (NonEmpty Name)] -> [Doc]
whereClause [] = []
whereClause as = [indent 2 $ "where" </> indent 2 (stack $ map defineLengthAttr as)]
  where
    defineLengthAttr a = renderHaskellAttrName a <+> "="
                            <+> "fromIntegral (length"
                            <+> renderHaskellName (NE.head $ attrInfo a)
                            <> ") :: Int64"

renderHaskellAttrName :: Attr a -> Doc
renderHaskellAttrName = renderHaskellName . attrName

functionBody :: ParsedOp -> Doc
functionBody pOp
    | parsedOpIsMonadic pOp
        = "build $ do"
            </> indent indentation (bindOpInputsVar
                        </> "buildOp" <+> outputListsSizes <+> opDef)
    | otherwise
        = "pureOp" <+> outputListsSizes <+> "$ do"
            </> indent indentation (bindOpInputsVar </> "return" <+> opDef)
  where
    outputListsSizes = brackets $ commasep
        [ renderHaskellName a
        | ParsedArg { parsedArgCase = ListArg { argLength = a } }
            <- parsedOutputs pOp
        ]
    opInputsVar = "op'inputs"
    bindOpInputsVar = opInputsVar <+> "<- fmap Prelude.concat $ Prelude.sequence"
                            <+> brackets (commasep $ map (\a -> "buildInputs" <+> a) tensorArgs)
    opDef = parens $ hang 0 $ stack $
        "opDef" <+> renderQuotedTFName (parsedOpName pOp) :
        -- Renders type parameter arguments.
        [ "& opAttr" <+> renderQuotedTFName n <+> ".~" <+> inferredTypeExpr a
        | a <- inferredTypeAttrs pOp, let n = attrName a
        ] ++
        -- Renders mandatory attributes as function parameters.
        [ "& opAttr" <+> renderQuotedTFName n <+> ".~" <+> renderHaskellName n
        | a <- explicitInputAttrs pOp, let n = attrName a
        ] ++
        -- Renders sizes of tensor list types having number_attr.
        [ "& opAttr" <+> renderQuotedTFName n <+> ".~" <+> renderHaskellName n
        | a <- inferredListSizeAttrs pOp, let n = attrName a
        ] ++
        ["& op'options & opInputs .~" <+> opInputsVar]
    tensorArgs = renderTensorArg <$> parsedInputs pOp
    renderTensorArg = renderHaskellName . parsedArgName
    inferredTypeExpr a
        | typeParamIsList $ attrInfo a
            = "fromTensorTypes (Proxy :: Proxy" <+> renderHaskellAttrName a
                    <> ")"
        | otherwise = "tensorType (undefined ::" <+> renderHaskellAttrName a
                            <> ")"

-- | Write a comment with the inputs/outputs/attributes in proto format, for
-- debugging.
extras :: OpDef -> Doc
extras d = enclose "{-\n" "\n-}" $
            strictText $ Text.pack $
            showMessage ((def :: OpDef)
                        & inputArg .~ (d ^. inputArg)
                        & outputArg .~ (d ^. outputArg)
                        & attr .~ (d ^. attr))

-- | The type signature for an op.
-- Of the form:
-- forall t1 t2 v1 v2 . (TensorType t1, TensorType t2)
--      => {pre} Float -> Tensor t1 v1 -> Tensor t2 v2
-- where "Float" is an explicit input attribute, "Tensor t1 v1" is an input, and
-- "Tensor t2 v2" is an output.
typeSig :: Doc -> ParsedOp -> Doc
typeSig pre pOp = constraints
            <+/> pre </> signatureFold (map attrInput (explicitInputAttrs pOp)
                                ++ map tensorArgAndComment (parsedInputs pOp)
                                ++ [outputs])
  where
    constraints
        | null classConstraints = empty
        | otherwise = "forall" <+> sep typeParams <+> "." <+> tuple classConstraints <+> "=>"
    typeParams = [strictText v | k <- parsedInputs pOp ++ parsedOutputs pOp,
                  ArgSomeTensor v <- [argKind $ parsedArgCase k]]
                ++ [renderHaskellAttrName n | n <- inferredTypeAttrs pOp]
                ++ if parsedOpIsMonadic pOp then ["m'"] else []
    -- Use m' as the type parameter to avoid clashing with an attribute name.
    monadConstraint
        | parsedOpIsMonadic pOp = ["MonadBuild m'"]
        | otherwise = []
    classConstraints = monadConstraint ++ map tensorArgConstraint
                                                    (inferredTypeAttrs pOp)
    signatureFold = folddoc (\x y -> x </> "->" <+> y)
    attrInput a = renderAttrType (attrInfo a) <+> hang 0 ("-- ^" <+> attrComment a)
    renderAttrType (AttrSingle a) = renderAttrBaseType a
    renderAttrType (AttrList a) = brackets $ renderAttrBaseType a
    renderAttrBaseType = \case
        AttrBytes -> "ByteString"
        AttrInt64 -> "Data.Int.Int64"
        AttrFloat -> "Float"
        AttrBool -> "Bool"
        AttrType -> "DataType"
        AttrShape -> "Shape"
        AttrTensor -> "TensorProto"

    tensorArgAndComment t = tensorArg t <+> hang 0 ("-- ^" <+> argComment t)
    outputs = case parsedOutputs pOp of
        [] -> wrapOutput "ControlNode"
        -- TODO(judahjacobson): To improve indentation: `tensorArgAndComment a`
        [a] -> wrapOutput (tensorArg a) <+> "-- ^" <+> argComment a
        as -> wrapOutput (tuple (map tensorArg as)) <+/> resultComment as
    wrapOutput o
        | parsedOpIsMonadic pOp = "m'" <+> parens o
        | otherwise = o

-- | Render an op input or output.
-- For example: "Tensor Ref Int64", "Tensor v t"
tensorArg :: ParsedArg -> Doc
tensorArg p = case parsedArgCase p of
    SimpleArg { argType = t, argKind = k } -> tensorType t k
    ListArg { argType = t, argKind = k } -> brackets $ tensorType t k
    MixedListArg {argTypeAttr = t, argKind = k}
        -> "TensorList" <+> parens (kind k) <+> renderHaskellName t
  where
    kind k = case k of
                ArgTensorRef -> "Ref"
                ArgTensorValue -> "Value"
                ArgTensorBuild -> "Build"
                ArgSomeTensor v -> strictText v
    tensorType t k = let
        a = case t of
                ArgTypeFixed dt -> strictText $ dtTypeToHaskell dt
                ArgTypeAttr n -> renderHaskellName n
        in "Tensor" <+> kind k <+> a

attrComment :: Attr a -> Doc
attrComment a = argComment' (attrName a) (attrDescription a)

argComment :: ParsedArg -> Doc
argComment a = argComment' (parsedArgName a) (parsedArgDescription a)

argComment' :: Name -> Text.Text -> Doc
argComment' argName argDesc =
    bold (renderTFName argName) <> splitMultilineText (":" <+>) argDesc

bold :: Doc -> Doc
bold n = "__" <> n <> "__"

-- | Comment for the outputs of an op.
-- For example:
--   -- ^ (__output1__, __output2__)
--   --
--   -- * __output1__: description1
--   --
--   -- * __output2__: description2
resultComment :: [ParsedArg] -> Doc
resultComment os = stack $ flatten commentSummary : map commentDetails os
  where
    commentSummary = "-- ^" <+> tuple [bold (renderTFName $ parsedArgName o) | o <- os]
    commentDetails o =
        stack [ "--"
              , "-- *" <+> argComment o
              ]

-- | Constraints for a given type parameter.
-- E.g.: "TensorType t" or "OneOf [Int64, Float] t"
-- or "TensorTypes ts" or "OneOfs [..] ts".
tensorArgConstraint :: Attr TypeParam -> Doc
tensorArgConstraint a = case attrInfo a of
    TypeParam False Nothing -> "TensorType" <+> n
    TypeParam False (Just as) -> "OneOf" <+> typeList as <+> n
    TypeParam True Nothing -> "TensorTypes" <+> n
    TypeParam True (Just as) -> "OneOfs" <+> typeList as <+> n
  where
    n = renderHaskellAttrName a
    -- Produces a type-level list, e.g.: '[Int32,Int64,Float]
    typeList = ("'" <>) . brackets . commasep . map strictText .
                    Set.toList . Set.fromList .
                    map dtTypeToHaskell . toList

-- NOTE: The cases of this function should be kept in sync with
-- TensorFlow.Types.AllTensorTypes.
dtTypeToHaskell :: DataType -> Text.Text
dtTypeToHaskell DT_BOOL = "Bool"
dtTypeToHaskell DT_BFLOAT16 = "Data.Word.Word16"
dtTypeToHaskell DT_COMPLEX128 = "(Data.Complex.Complex Double)"
dtTypeToHaskell DT_COMPLEX64 = "(Data.Complex.Complex Float)"
dtTypeToHaskell DT_DOUBLE = "Double"
dtTypeToHaskell DT_FLOAT = "Float"
dtTypeToHaskell DT_INT16 = "Data.Int.Int16"
dtTypeToHaskell DT_INT32 = "Data.Int.Int32"
dtTypeToHaskell DT_INT64 = "Data.Int.Int64"
dtTypeToHaskell DT_INT8 = "Data.Int.Int8"
dtTypeToHaskell DT_QINT32 = "Data.Int.Int32"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_QINT8 = "Data.Word.Word8"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_QINT16 = "Data.Int.Int16"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_QUINT16 = "Data.Word.Word16"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_QUINT8 = "Data.Word.Word8"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_STRING = "Data.ByteString.ByteString"
dtTypeToHaskell DT_UINT16 = "Data.Word.Word16"
dtTypeToHaskell DT_UINT32 = "Data.Word.Word32"
dtTypeToHaskell DT_UINT64 = "Data.Word.Word64"
dtTypeToHaskell DT_HALF = "Data.Word.Word16"  -- TODO(gnezdo): make unique
dtTypeToHaskell DT_UINT8 = "Data.Word.Word8"
dtTypeToHaskell DT_RESOURCE = "ResourceHandle"
dtTypeToHaskell DT_VARIANT = "Variant"
dtTypeToHaskell x =
    Text.pack $ "Unsupported type in dtTypeToHaskell: " ++ show x

-- | haddockComment escapes TensorFlow doc strings into haddock.
-- TODO(gnezdo): deal with the markup.
haddockComment :: Text.Text -> Doc
haddockComment = strictText

-- | Generate a multiline comment.  For example:
--   summary'
--   --
--   -- detail_line1
--   -- detail_line2
--   -- ...
multilineComment :: Text.Text -> Text.Text -> Doc
multilineComment summary' detail =
    haddockComment summary' </>
    splitMultilineText insertParagraphAndComment detail
  where insertParagraphAndComment x = "--" </> "--" <+> x

-- | Converts the given multi-line detail string into
-- a multi-line haddock. Applies the given lead to the
-- first line. Returns an empty document for empty detail.
splitMultilineText :: (Doc -> Doc) -> Text.Text -> Doc
splitMultilineText lead detail =
  case Text.lines detail of
    [] -> empty
    (l : ls) -> stack $ lead (haddockComment l)
                      : map (("--" <+>) . haddockComment) ls

indentation :: Int
indentation = 4


text:
module NeuralNetwork
  () where

import           Control.Monad         (mapM)
import qualified Data.Sequence         as Seq
import           Debug.Trace
import           Numeric.LinearAlgebra ((<>), (><))
import qualified Numeric.LinearAlgebra as LA

type MatrixNum = LA.Matrix Double
type Cache = ((MatrixNum, MatrixNum, MatrixNum), MatrixNum)

data Parameter = Parameter
  { weightsParam :: MatrixNum
  , biasesParam  :: MatrixNum
  }

-- Utils
emptyx :: MatrixNum
emptyx = (0><0) []

emptyCache :: Cache
emptyCache = ((emptyx, emptyx, emptyx), emptyx)

-- Activation functions
relu :: MatrixNum -> MatrixNum
relu = LA.cmap (max 0)

sigmoid :: MatrixNum -> MatrixNum
sigmoid = LA.cmap (\x -> (1 / (1 + exp (-1 * x))))

relu' :: MatrixNum -> MatrixNum
relu' = LA.cmap f
  where
    f x | x < 0 = 0
        | otherwise = 1

sigmoid' :: MatrixNum -> MatrixNum
sigmoid' xw = (sigmoid xw) * (1 - sigmoid xw)

-- | Initialize parameters. Randomize weights and zeroes biases
initParams :: [Int] -> IO [Parameter]
initParams layerDims =
  (fmap . fmap) (\(bx, wx) -> Parameter wx bx) ((zip biases) <$> weights)
    where
      weights = mapM genParams layerDims'
      layerDims' = (zip layerDims (tail layerDims))
      genParams (prevDim, currentDim) = LA.randn currentDim prevDim
      biases =
        fmap
        (\(_, currentDim) -> (currentDim><1) (replicate currentDim 0))
        layerDims'

-- | Propagate forward to next layer. Calulates next a
linearForward axPrev wx bx activationFun =
  let zx = linearForward' axPrev wx bx
      ax = activationFun zx
      cache = ((axPrev, wx, bx), zx)
      linearForward' axPrev wx bx = (wx <> axPrev) + bx :: MatrixNum
  in  (ax, cache)

-- | Complete forward propagation. Calculates final ax
propagateForward :: MatrixNum
                 -> [Parameter]
                 -> (MatrixNum, [Cache])
propagateForward xx params  =
  (outputAx, hiddenCache ++ [outputCache])
    where
      params' = take ((length params) - 1) params
      f (aPrev, _) (Parameter wx bx) = linearForward aPrev wx bx relu
      hidden = scanl f (xx, emptyCache) params'
      hiddenCache = drop 1 (fmap snd hidden)
      (lastHiddenAx, _) = last hidden
      Parameter wxLast bxLast = last params
      (outputAx, outputCache) = linearForward lastHiddenAx wxLast bxLast sigmoid

-- | Propagate backward to previous layer
linearBackward :: MatrixNum -> Cache -> (MatrixNum -> MatrixNum)
               -> (MatrixNum, MatrixNum, MatrixNum)
linearBackward dax ((axPrev, wx, bx), zx) af' =
  (daxPrev, dwx, dbx)
    where
      dzx = dax * af' zx
      m = fromIntegral $ (snd . LA.size) axPrev
      dwx = (dzx <> (LA.tr axPrev)) / m
      sums = fmap sum (LA.toLists dzx)
      layerSize = (fst . LA.size) wx
      dbx = ((layerSize><1) sums) / m
      daxPrev = (LA.tr wx) <> dzx

-- | Complete backpropagation. Computes gradients (dax, dwx, dbx)
-- of each layer
propagateBackward :: MatrixNum -> (LA.Vector Double) -> [Cache]
                  -> [(MatrixNum, MatrixNum, MatrixNum)]
propagateBackward axOut yx caches =
  let l = length caches
      yx' = LA.reshape 1 yx
      daxOut = - (yx' / axOut) - ((1 - yx') / (1 - axOut))
      gradOut = linearBackward daxOut (last caches) sigmoid'
      hiddenCaches = reverse $ take (length caches - 1) caches
      grads = scanl f gradOut hiddenCaches
      f (dax, _, _) cache = linearBackward dax cache relu'
  in  reverse grads

-- | Perform gradient descent. Computes new decreased gradients (dax, dwx, dbx)
updateParameters :: [Parameter]
                 -> [(MatrixNum, MatrixNum, MatrixNum)]
                 -> Double
                 -> [Parameter]
updateParameters params grads eta = fmap fun (zip params grads)
   where fun ((Parameter wx bx), (_, dwx, dbx)) =
           Parameter
             (wx - (LA.scalar eta) * dwx)
             (bx - (LA.scalar eta) * dbx)

train :: MatrixNum
      -> LA.Vector Double
      -> [Int]
      -> Double
      -> Int
      -> IO [Parameter]
train xx yx layerDims eta iterations =
  do
    params <- initParams layerDims
    return $ train' params iterations
  where
    train' params 0 = params
    train' params times =
      let (ax, caches) = propagateForward xx params
          grads = propagateBackward ax yx caches
          newParams = updateParameters params grads eta
      in  train' newParams (times - 1)

max_stars_repo_name: abs-zero/haskellml

max_issues_repo_name: abs-zero/haskellml

max_forks_repo_name: abs-zero/haskellml


text:
{-# LANGUAGE RecordWildCards #-}

import Control.Monad.Trans.Except
import Data.Word
import Data.Complex
import Data.Maybe
import Data.Monoid

import Control.Error.Util
import Pipes                as P
import Pipes.Prelude        as P
import Options.Applicative
import Data.Vector.Storable as VS hiding ((++))
import Data.Vector.Generic  as VG hiding ((++))

import SDR.Util             as U
import SDR.RTLSDRStream
import SDR.Filter
import SDR.ArgUtils
import SDR.FilterDesign
import SDR.Pulse
import SDR.PipeUtils
import SDR.CPUID

data Options = Options {
    frequency    :: Word32,
    bandwidth    :: Maybe Int
}

optParser :: Parser Options
optParser = Options
          <$> option (fmap fromIntegral parseSize) (
                 long "frequency"
              <> short 'f'
              <> metavar "FREQUENCY"
              <> help "Frequency to tune to"
              )
          <*> optional (option (fmap fromIntegral parseSize) (
                 long "bandwidth"
              <> short 'b'
              <> metavar "BANDWIDTH"
              <> help "Filter bandwidth. From 0 to 32K."
              ))

opt :: ParserInfo Options
opt = info (helper <*> optParser) (fullDesc <> progDesc "Receive AM Radio" <> header "RTLSDR AM")

size = 4096

doIt Options{..} = do

    info            <- lift getCPUInfo

    str             <- sdrStream (defaultRTLSDRParams (frequency + 256000) 1024000) 1 (fromIntegral $ size * 2)

    let coeffsDecim :: [Float]
        coeffsDecim =  VS.toList $ windowedSinc 71 0.4 blackman
    deci            <- lift $ fastDecimatorC info 2 coeffsDecim

    let coeffsFilt  :: [Float]
        coeffsFilt  =  VS.toList $ windowedSinc 71 (fromIntegral (fromMaybe 16000 bandwidth) / 32000) blackman
    filt            <- lift $ fastFilterC info coeffsFilt

    let coeffsResp  :: [Float]
        coeffsResp  =  VS.toList $ windowedSinc 71 0.25 blackman
    resp            <- lift $ fastResamplerR info 3 2 coeffsResp

    pulseSink       <- lift pulseAudioSink

    lift $ runEffect $   str
                     >-> P.map (interleavedIQUnsignedByteToFloatFast info)
                     >-> P.map (VG.zipWith (*) (quarterBandUp size))
                     >-> firDecimator deci size
                     >-> firDecimator deci size
                     >-> firDecimator deci size
                     >-> firDecimator deci size
                     >-> firDecimator deci size
                     >-> firFilter  filt size
                     >-> P.map (VG.map magnitude)
                     >-> dcBlockingFilter
                     >-> firResampler resp size
                     >-> P.map (VG.map (* 6))
                     >-> pulseSink

main = execParser opt >>= exceptT putStrLn return . doIt

max_stars_repo_name: adamwalker/sdr-apps

max_issues_repo_name: adamwalker/sdr-apps

max_forks_repo_name: adamwalker/sdr-apps


text:
module HVX.Internal.Matrix
  ( Mat
  , allMat
  , anyMat
  , diagMat
  , ei
  , fpequalsMat
  , lpnorm
  , matrixPow
  , reduceMat
  , scalarMat
  , zeroMat
  , zeroVec
  ) where

import Numeric.LinearAlgebra hiding (i)
import Numeric.LinearAlgebra.Util

import HVX.Internal.Util

type Mat = Matrix Double

allMat :: (Double -> Bool) -> Mat -> Bool
allMat f x = all f (toList . flatten $ x)

anyMat :: (Double -> Bool) -> Mat -> Bool
anyMat f x = any f (toList . flatten $ x)

diagMat :: Mat -> Mat
diagMat = diag . flatten

ei :: Int -> Int -> Mat
ei n i = buildMatrix n 1 (\(j, _) -> if i == j then 1 else 0)

fpequalsMat :: Mat -> Mat -> Bool
fpequalsMat a b
  | ra == rb && ca == cb = all (uncurry fpequals) $ zip alist blist
  | otherwise = error "Two matrices with different dimensions cannot possibley be equal!"
  where
    ra = rows a
    rb = rows b
    ca = cols a
    cb = cols b
    alist = toList . flatten $ a
    blist = toList . flatten $ b

lpnorm :: Double -> Mat -> Double
lpnorm p x = sumElements y ** (1/p)
  where pMat = (1><1) [p]
        y = abs x ** pMat

matrixPow :: Double -> Mat -> Mat
matrixPow p = mapMatrix (** p)

reduceMat :: ([Double] -> Double) -> Mat -> Mat
reduceMat f = (1><1) . (:[]) . f . toList . flatten

scalarMat :: Double -> Mat
scalarMat x = (1><1) [x]

zeroMat :: Int -> Mat
zeroMat n = zeros n n

zeroVec :: Int -> Mat
zeroVec n = zeros n 1

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
{-# LANGUAGE ExistentialQuantification , DeriveFunctor,FlexibleContexts,ImpredicativeTypes,TypeFamilies #-}

module Exponential.SO3 where

import Exponential.Class
import Data.Complex

import Control.Lens

import Linear.Vector
import Linear.V3
import Linear.Affine
import Linear.Matrix
import Linear.Metric
import Local
import Linear.V4

import Rotation.SO3
import Space.Class
import qualified Debug.Trace as T

instance Exponential SO3 where
    logM = logM3 . unSO3
    expM = SO3 . expM3

data CSO3 a = CSO3 { unCSO3 :: V3 (V3 a)} deriving(Show,Read,Functor)

mkLorentz :: Num a => M33 a -> V3 a -> M44 a
mkLorentz (V3 r1 r2 r3) (V3 tx ty tz) =
  V4 (snoc3 r1 tx) (snoc3 r2 ty) (snoc3 r3 tz) (V4 tx ty tz 0)
  where snoc3 (V3 x y z) = V4 x y z

expM31 v1@(V3 ox oy oz) v2@(V3 dx dy dz)  = foldr1 (!+!) $ zipWith (*!!) (fmap realPart [l1,l2,l3,l4]) [eye4,a , a!*!a , a!*!a !*!a]
    where
        alpha = v1 `dot` v1  - v2 `dot` v2
        beta = v1 `dot` v2
        delta = sqrt $  alpha^2 + 4*beta
        lambda1 = 0 :+ sqrt (abs((alpha - delta )/2))
        lambda2 =  sqrt ((alpha + delta )/2) :+  0
        l1 = (-cosh(lambda1)*lambda2^2 + lambda1^2*cosh(lambda2))/(lambda1^2 - lambda2^2)
        l2 = (-sinh(lambda1)*lambda2^3 + lambda1^3*sinh(lambda2))/(lambda1^2 - lambda2^2)/lambda1/lambda2
        l3 = (cosh lambda1  + cosh lambda2 )/(lambda1^2 - lambda2^2)
        l4 = (lambda2^2 * sinh lambda1 - lambda1* sinh lambda2)/(lambda1^2 - lambda2^2)/lambda1/lambda2
        a = mkLorentz (skewV3 v1) v2

expM3 v@(V3 x y z)
    = V3
        (V3 (cos t + c*x^2) (-s*z +c*x*y) (s*y+c*x*z))
        (V3 (s*z+c*x*y) (cos t + c*y^2) (-s*x+c*y*z))
        (V3 (-s*y+c*x*z ) (s*x+c*y*z) (cos t + c*z^2))
    where
        t = norm v
        c =  (1 - cos t)/(t^2 + 1.24e-16)
        s = sinc t

logM3 :: (Floating a, Ord a) => (V3 (V3 a)) -> V3 a
logM3 v = (t /( 2*sin t + 1.24e-16)) *^ (V3 (v^._z._y - v^._y._z) (v^._x._z - v^._z._x) (v^._y._x - v^._x._y))
    where
        t = acosc $ (trace v - 1)/2

acosc x
 | x > 1 = acos 1
 | x < -1 = acos(-1)
 | otherwise = acos x

sinc x
 | x == 0 = 1
 | otherwise = sin x / x

max_stars_repo_name: massudaw/mtk

max_issues_repo_name: massudaw/mtk

max_forks_repo_name: massudaw/mtk


text:
module Tests.Parametric (tests) where

import Data.Maybe (fromJust)
import Statistics.Test.StudentT
import Statistics.Types
import qualified Data.Vector.Unboxed as U
import Test.Tasty (testGroup)
import Tests.Helpers  (testEquality)
import qualified Test.Tasty as Tst

tests :: Tst.TestTree
tests = testGroup "Parametric tests" studentTTests

-- 2 samples x 20 obs data
--
-- Both samples are samples from normal distributions with the same variance (= 1.0),
-- but their means are different (0.0 and 0.5, respectively).
--
-- You can reproduce the data with R (3.1.0) as follows:
--   set.seed(0)
--   sample1 = rnorm(20)
--   sample2 = rnorm(20, 0.5)
--   student = t.test(sample1, sample2, var.equal=T)
--   welch = t.test(sample1, sample2)
--   paired = t.test(sample1, sample2, paired=T)
sample1, sample2 :: U.Vector Double
sample1 = U.fromList [
  1.262954284880793e+00,
 -3.262333607056494e-01,
  1.329799262922501e+00,
  1.272429321429405e+00,
  4.146414344564082e-01,
 -1.539950041903710e+00,
 -9.285670347135381e-01,
 -2.947204467905602e-01,
 -5.767172747536955e-03,
  2.404653388857951e+00,
  7.635934611404596e-01,
 -7.990092489893682e-01,
 -1.147657009236351e+00,
 -2.894615736882233e-01,
 -2.992151178973161e-01,
 -4.115108327950670e-01,
  2.522234481561323e-01,
 -8.919211272845686e-01,
  4.356832993557186e-01,
 -1.237538421929958e+00]
sample2 = U.fromList [
  2.757321147216907e-01,
  8.773956459817011e-01,
  6.333363608148415e-01,
  1.304189509744908e+00,
  4.428932256161913e-01,
  1.003607972233726e+00,
  1.585769362145687e+00,
 -1.909538396968303e-01,
 -7.845993538721883e-01,
  5.467261721883520e-01,
  2.642934435604988e-01,
 -4.288825501025439e-02,
  6.668968254321778e-02,
 -1.494716467962331e-01,
  1.226750747385451e+00,
  1.651911754087200e+00,
  1.492160365445798e+00,
  7.048689050811874e-02,
  1.738304100853380e+00,
  2.206537181457307e-01]

testTTest :: String
          -> PValue Double
          -> Test d
          -> [Tst.TestTree]
testTTest name pVal test =
  [ testEquality name (isSignificant pVal test) NotSignificant
  , testEquality name (isSignificant (mkPValue $ pValue pVal + 1e-5) test)
    Significant
  ]

studentTTests :: [Tst.TestTree]
studentTTests = concat
  [ -- R: t.test(sample1, sample2, alt="two.sided", var.equal=T)
    testTTest "two-sample t-test SamplesDiffer Student"
      (mkPValue 0.03410) (fromJust $ studentTTest SamplesDiffer sample1 sample2)
    -- R: t.test(sample1, sample2, alt="two.sided", var.equal=F)
  , testTTest "two-sample t-test SamplesDiffer Welch"
      (mkPValue 0.03483) (fromJust $ welchTTest SamplesDiffer sample1 sample2)
    -- R: t.test(sample1, sample2, alt="two.sided", paired=T)
  , testTTest "two-sample t-test SamplesDiffer Paired"
      (mkPValue 0.03411) (fromJust $ pairedTTest SamplesDiffer sample12)
    -- R: t.test(sample1, sample2, alt="less", var.equal=T)
  , testTTest "two-sample t-test BGreater Student"
      (mkPValue 0.01705) (fromJust $ studentTTest BGreater sample1 sample2)
    -- R: t.test(sample1, sample2, alt="less", var.equal=F)
  , testTTest "two-sample t-test BGreater Welch"
      (mkPValue 0.01741) (fromJust $ welchTTest BGreater sample1 sample2)
    -- R: t.test(sample1, sample2, alt="less", paired=F)
  , testTTest "two-sample t-test BGreater Paired"
      (mkPValue 0.01705) (fromJust $ pairedTTest BGreater sample12)
  ]
  where sample12 = U.zip sample1 sample2


text:
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE TypeFamilies #-}

-- Import all models under maintenance.
-- Models not imported here will not be compiled
-- when invoking `stack bench`.
import qualified BetaBin
import Control.Applicative ((<$>))
import Control.Arrow (first, second)
import Control.Monad (unless, when)
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Enumerator
import Control.Monad.Bayes.Inference
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Simple
import Control.Monad.Bayes.Weighted
import Control.Monad.IO.Class (MonadIO, liftIO)
import qualified DPmixture
import qualified Data.Vector as Vector
import qualified Dice
import qualified Gamma
import Graphics.Rendering.Chart.Backend.Cairo
import Graphics.Rendering.Chart.Easy
import qualified HMM
import Numeric.LogDomain
import Options.Applicative
import Plotting
import Statistics.Sample
import System.Directory
import System.IO

opts :: ParserInfo Bool
opts =
  flip info fullDesc $
    switch
      ( long "trial"
          <> help "Run a quick version of benchmarks to check that all is working correctly"
      )

tryCache :: (MonadIO m, Read a, Show a) => FilePath -> m a -> m a
tryCache filepath fresh = do
  exists <- liftIO $ doesFileExist filepath
  if exists
    then liftIO $ fmap read (readFile filepath)
    else do
      value <- fresh
      liftIO $ writeFile filepath (show value)
      return value

main = do
  -- make sure `putStrLn` prints to console immediately
  hSetBuffering stdout LineBuffering
  trial <- execParser opts
  when trial $ putStrLn "Trial run"
  sampleIO hmmBenchmark

meanVar :: (MonadDist m, CustomReal m ~ Double) => Int -> m Double -> m (Double, Double)
meanVar n d = do
  xs <- Vector.replicateM n d
  return $ meanVariance xs

smcParams :: [Int]
smcParams = [10, 20, 50, 100, 200, 500, 1000]

ns :: [Int]
--ns = [10,20,50,100,200,500,1000]
ns = [10, 20 .. 1000]

smcParamsDouble :: [Double]
smcParamsDouble = map fromIntegral smcParams

smcResults :: (MonadDist m, CustomReal m ~ Double) => [m (Vector.Vector Double)]
smcResults = map (\p -> Vector.replicateM 10 $ fmap HMM.hmmKL $ explicitPopulation $ smcMultinomial (length HMM.values) p HMM.hmm) smcParams

hmmBenchmark :: SamplerIO ()
hmmBenchmark = do
  liftIO $ putStrLn "running HMM benchmark"
  isSamples <- fmap (drop 5000) $ explicitPopulation $ importance 10000 HMM.hmm
  let isRes = map (\n -> HMM.hmmKL $ take n isSamples) ns
  -- mhSamples <- fmap (drop 5000) $ traceMH 10000 HMM.hmm
  -- let mhRes = map (\n -> HMM.hmmKL $ take n $ map (,1) mhSamples) ns
  mhPriorSamples <- drop 5000 <$> mhPrior HMM.hmm 10000
  let mhPriorRes = map (\n -> HMM.hmmKL $ take n $ map (,1) mhPriorSamples) ns
  pimhSamples <- pimh (length HMM.values) 100 1000 HMM.hmm
  let pimhRes = map (\n -> HMM.hmmKL $ take n $ map (,1) pimhSamples) ns
  liftIO $ toFile (fo_format .~ PDF $ def) "anytime.pdf" $ do
    layout_title .= "HMM"
    anytimePlot
      "#samples"
      "KL"
      ns
      [ ("IS", isRes),
        -- ("MHtrace", mhRes),
        ("MHprior", mhPriorRes),
        ("PIMH", pimhRes)
      ]
  smcRes <- sequence smcResults
  liftIO $ toFile (fo_format .~ PDF $ def) "smc.pdf" $ do
    layout_title .= "HMM"
    oneShotPlot
      "#particles"
      "KL"
      [ ("SMC", zip smcParamsDouble (map (errBars 2) smcRes))
      ]

max_stars_repo_name: saeedhadikhanloo/monad-bayes

max_issues_repo_name: saeedhadikhanloo/monad-bayes

max_forks_repo_name: saeedhadikhanloo/monad-bayes


text:
module MachineLearning.NeuralNetwork.TopologyTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.NeuralNetwork.Topology
import qualified MachineLearning.NeuralNetwork.TopologyMaker as TM

nnt = TM.makeTopology TM.ASigmoid TM.LLogistic 15 2 [10]

flattenTest = do
  theta <- initializeThetaIO nnt
  let theta' = flatten $ unflatten nnt theta
      norm = LA.norm_2 (theta - theta')
  assertApproxEqual "flatten" 1e-10 0 norm

tests = [ testGroup "flatten" [
            testCase "flatten" flattenTest
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeSynonymInstances #-}
module Statistics.BBVI.Distribution.Dirichlet
  ( alphas
  , dirichlet
  , Dirichlet(..)
  )
where

import qualified Data.Vector                   as V
import           Numeric.SpecFunctions          ( logGamma
                                                , digamma
                                                )
import           Statistics.BBVI.Class
import           Statistics.BBVI.Propagator     ( SampleVector )
import qualified System.Random.MWC.Distributions
                                               as MWCD

-- | data type for a dirichlet distribution
newtype Dirichlet = Diri (V.Vector Double) deriving (Show, Eq, Ord, Read)

-- | build a dirichlet
dirichlet
  :: V.Vector Double -- ^ vector of parameters
  -> Dirichlet
dirichlet xs = Diri $ V.map (max 1e-10) xs

-- | get dirichlet parameters
alphas :: Dirichlet -> V.Vector Double -- ^ dirichlet parameters
alphas (Diri xs) = xs

logB :: Dirichlet -> Double
logB d = V.sum (V.map logGamma as) - logGamma (V.sum as) where as = alphas d

instance DistUtil Dirichlet where
  nParams         = V.length . alphas
  fromParamVector = dirichlet
  toParamVector   = alphas

instance Dist Dirichlet SampleVector where
  resample d gen = MWCD.dirichlet (alphas d) gen
  logProb d cat = V.sum (V.zipWith (\alpha x -> (alpha - 1) * log x) as cat)
    - logB d
    where as = alphas d
  paramGradOfLogQ d cat = V.zipWith (\a x -> summed - digamma a + log x) as cat
   where
    as     = alphas d
    summed = digamma (V.sum as)

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
{-# OPTIONS -O3 -optc-O3 -XFlexibleInstances #-}

module Paraiso(
  Architecture(..),
  Dsl,
  Expr(..),
  Register(..),(=$),
  compile,
  parallel,
  sequential,
  cuda,
  mesh
)where

import Control.Monad.RWS
import Data.Complex
import Data.List
import Data.Maybe

import Util
import Debug

type VarID = String
type VarModifier = String -> String
modRegister::VarModifier
modRegister = (++ "_reg")
modDevice::VarModifier
modDevice = (++ "_dev")
flagRegister = ""  --"__shared__ "

data Architecture = X86 |
                    CUDA{numberOfGrids::Int, numberOfThreads::Int} |
                    MPI Architecture deriving Eq

data Statement = Statement{funcName::String, phase::StatementType}

data StatementType = Alloc String | Store String | Free String |
              FuncCall String [String] [String]
              deriving (Eq, Ord, Show, Read)
funcNameMain = "main"

data HardwareContext = Host | Device String deriving (Eq, Ord, Show, Read)
data Context = Context
    {
     precision::String,
     freeVariable::VarID, storeModifier::VarModifier, allocModifier::VarModifier,
     parallelDegree::Int,
     architecture:: Architecture,
     hardware::HardwareContext,
     fogArgument::[String],
     fogArgumentType::[String],
     openHostParallel::[String],
     closeHostParallel::[String],
     memSendHD::[String],
     memRecvHD::[String],
     memSendDH::[String],
     memRecvDH::[String]
    }

data Expr a =
    Imm a |
    Rand a a |
    Var VarID |
    Add (Expr a) (Expr a) |
    Sub (Expr a) (Expr a) |
    Mul (Expr a) (Expr a) |
    Div (Expr a) (Expr a) |
    Abs (Expr a) |
    Rec (Expr a) |
    Sgn (Expr a) deriving (Eq, Ord, Show, Read)

compileExpr:: (Show a) => VarModifier -> Expr a -> String
compileExpr mo = c where
    c (Imm a) = show a
    c (Rand lo hi) = "drand(" ++ show lo ++ " , " ++ show hi ++ ")"
    c (Var v) = mo v
    c (Add a b) = "(" ++ c a ++ " + " ++ c b ++ ")"
    c (Sub a b) = "(" ++ c a ++ " - " ++ c b ++ ")"
    c (Mul a b) = "(" ++ c a ++ " * " ++ c b ++ ")"
    c (Div a b) = "(" ++ c a ++ " / " ++ c b ++ ")"
    c (Abs a) = "abs(" ++ c a ++ ")"
    c (Rec a) = "(1.0 / " ++ c a ++ ")"
    c (Sgn a) = "(" ++ c a ++ " >0?1:-1)"

instance (Num t) => Num (Expr t) where
    a + b = Add a b
    a - b = Sub a b
    a * b = Mul a b
    abs a = Abs a
    signum a = Sgn a
    fromInteger i =  Imm (fromInteger i)

instance Fractional t => Fractional (Expr t) where
    a / b = Div a b
    recip a = Rec a
    fromRational i = Imm (fromRational i)

instance Real t => Real (Expr t) where
    toRational = undefined

instance Floating t => Floating (Expr t) where
    pi = undefined
    exp = undefined
    log = undefined
    sin = undefined
    cos = undefined
    asin = undefined
    atan = undefined
    acos = undefined
    sinh = undefined
    cosh = undefined
    asinh = undefined
    atanh = undefined
    acosh = undefined

instance RealFloat t => RealFloat (Expr t) where
    floatRadix = undefined
    floatDigits = undefined
    floatRange = undefined
    decodeFloat = undefined
    encodeFloat = undefined
    isNaN = undefined
    isInfinite = undefined
    isDenormalized = undefined
    isNegativeZero = undefined
    isIEEE = undefined

instance RealFrac t => RealFrac (Expr t) where
    properFraction = undefined

type Dsl a = RWS () [Statement] Context a

tellA::[String] -> Dsl ()
tellA xs = do
  s <- get
  case hardware s of
    Host -> tellHA xs
    Device fn -> tellDA fn xs
tellS::[String] -> Dsl ()
tellS xs = do
  s <- get
  case hardware s of
    Host -> tellHS xs
    Device fn -> tellDS fn xs

tellHA :: [String] -> Dsl ()
tellHA = tell.map (\st -> Statement funcNameMain (Alloc st))
tellHS :: [String] -> Dsl ()
tellHS = tell.map (\st -> Statement funcNameMain (Store st))
tellHF :: [String] -> Dsl ()
tellHF = tell.map (\st -> Statement funcNameMain (Free st))
tellDA :: String -> [String] -> Dsl ()
tellDA fn = tell.map (\st -> Statement fn (Free st))
tellDS :: String -> [String] -> Dsl ()
tellDS fn = tell.map (\st -> Statement fn (Free st))

getVarName = do
  let inc (h:ts) = (h:(show.(1+).read $ ts))
  state <- get
  let v = freeVariable state
  put state{freeVariable = inc v}
  return v

class Register a where
    allocate :: Dsl a
    store::a -> a -> Dsl ()
    output::[a] -> Dsl ()

(=$) :: Register a => a -> a -> Dsl ()
(=$) = store
infix 0 =$

instance Register (Expr Double) where
    allocate = do
      v <- getVarName
      s <- get
      let
        am = allocModifier s
        sm = storeModifier s
        prec = precision s
        arch = architecture s
        parD = parallelDegree s
      if hardware s == Host
         then do
           case arch of
             X86 -> tellA [am v  ++";"]
             CUDA{} -> do
               let
                 vDev = modDevice v
                 vSha  = modRegister v
                 sos = "sizeof(" ++ prec ++ ")*" ++ show parD
               tellHA  [prec ++ " *" ++ v ++ " = (" ++prec ++ "*) malloc(" ++sos++ ");"]
               tellHA  [prec ++ " *" ++ vDev ++ ";"]
               tellHA  ["cudaMalloc((void**) &" ++ vDev ++" ," ++ sos ++ ");"]
               s<-get
               put s{
                 memSendHD = (memSendHD s +:) $
                   "cudaMemcpy(" ++ vDev ++ " , " ++ v ++ ","
                      ++ sos ++ " , cudaMemcpyHostToDevice);",
                 memRecvHD = (memRecvHD s ++) $
                   [prec ++ " "++ flagRegister ++" " ++ vSha ++ ";" , vSha ++ " = " ++ sm vDev ++ ";"],
                 memSendDH = (memSendDH s +:) $
                   sm vDev ++ " = " ++ vSha ++ ";",
                 memRecvDH = (memRecvDH s +:) $
                   "cudaMemcpy(" ++ v ++ " , " ++ vDev ++ ","
                     ++ sos ++ ", cudaMemcpyDeviceToHost);",
                 fogArgument     = (fogArgument     s +:) $ vDev,
                 fogArgumentType = (fogArgumentType s +:) $ prec ++  " *" ++  vDev
               }
         else do
           tellA [prec ++ " " ++  am v ++";"]
      return $ Var v
    store addr expr = do
      case addr of
        Var v -> do
                  s <- get
                  let sm = storeModifier s
                  tellS [sm v ++ " = " ++ compileExpr sm expr ++ ";"]
        _ -> error "lhs is not a variable address."
    output addrs = do
      s <- get
      let
          sm = storeModifier s
          fromVar ex = case ex of
                       Var id -> sm id
                       _      -> "output is not a variable address."
          statement =
              unwords $
              (++ ["<< endl;"] ) $
              (["cout <<"] ++ ) $
              intersperse ("<< \" \" <<") $
              map fromVar addrs
      tellS [statement]

instance (RealFloat a, Register a) => Register (Complex a) where
    allocate = do
      r <- allocate
      i <- allocate
      return (r:+i)
    store (tr:+ti) (sr:+si) = do
      tmp <- allocate
      store tmp sr
      store ti si
      store tr tmp
    output = undefined

compile :: Architecture -> Dsl a-> String
compile arch a = code
  where
      (_,s,w) = runRWS a () initContext

      initContext = Context
            {
             freeVariable = "a1",
             storeModifier = id,
             allocModifier = ((precision initContext ++ " ") ++),
             parallelDegree = 1,
             architecture = arch,
             hardware = Host,
             precision = case arch of
                           CUDA{} -> "float"
                           _      -> "double",
             fogArgument =[],
             fogArgumentType =[],
             openHostParallel = [],
             closeHostParallel = [],
             memSendHD = [],
             memRecvHD = [],
             memSendDH = [],
             memRecvDH = []
            }

      funcAndBodies = groupSort funcName $ w
      funcCalls::[StatementType]
      funcCalls = concatMap (\st -> case st of
                                   (Statement _ fc@FuncCall{}) -> [fc]
                                   _ -> []) w
      lookupFunc fn = case find (\(FuncCall n _ _) -> n==fn) funcCalls of
                 Just fc -> fc
                 Nothing -> error $ "function disappeared! " ++ fn

      code = unlines $ progHeader ++ concatMap buildFunc funcAndBodies

      buildFunc (fn, stmts) = let
          (header,footer) = sandwich fn
          mysort = sortBy (\a b->f a `compare` f b)
          f (Statement _ Alloc{})    = 1
          f (Statement _ Store{})    = 2
          f (Statement _ FuncCall{}) = 2
          f (Statement _ Free{})     = 3
        in header ++ (map peal $ mysort stmts) ++ footer

      peal (Statement _ (Alloc s)) = s
      peal (Statement _ (Store s)) = s
      peal (Statement _ (Free  s)) = s
      peal (Statement _ (FuncCall fn args argtypes))
          = fn ++ "<<<grids,threads>>>(" ++ (unwords . intersperse ",") args ++ ");"

      progHeader = includes ++ utilfuncs
      includes =
        (++["using namespace std;"]) $
        map (\fn -> "#include <" ++ fn ++ ">") $
        headers
      headers = ["iostream", "cstdlib"]
        ++ case arch of
             CUDA{} -> [] -- ["cutil.h"]
             _      -> []
      utilfuncs =
        [
         "double drand(double lo, double hi){",
         "return lo + rand()/(double)RAND_MAX * (hi-lo);",
         "}"
        ]

      sandwich fn = if fn == funcNameMain then
           (["int main(int argc, char **argv){"]
            ++ case arch of
              CUDA{numberOfGrids = ng,  numberOfThreads = nt}
                -> ["dim3 grids(" ++ show ng ++");",
                         "dim3 threads(" ++ show nt ++");"]
              _ -> []

           , ["return 0;","}"])
        else
            (["__global__ void " ++ fn ++ "(" ++ (unwords $ intersperse "," $ att fn) ++ "){"]
            , ["}"])

      --argumentTypeTable
      att fn = let (FuncCall _ _ fat) = lookupFunc fn in fat

parallel size dsl = do
  iv <- getVarName
  state <- get
  put $ state{
              storeModifier = (++ ("[" ++ iv ++ "]")) . storeModifier state,
              allocModifier = (++ ("[" ++ show size ++ "]")) . allocModifier state,
              memRecvHD = (:memRecvHD state) $ "int " ++ iv ++ " = blockIdx.x * gridDim.x + threadIdx.x;"
--              fogArgument = (:fogArgument state) $  iv, NEEDED_FOR_GREATER_LOOP
--              fogArgumentType = (:fogArgumentType state) $ "int " ++ iv
             }
  let
    open  = "for(int "++iv++" = 0 ; "++iv++" < " ++ show size ++ " ; ++"++iv++"){"
    close = "}"
  tellS [open]
  do
    state <- get
    put $ state{
            openHostParallel  = openHostParallel state ++ [open],
            closeHostParallel = [close] ++ closeHostParallel state,
            parallelDegree = size * parallelDegree state
          }
    dsl
  tellS [close]
  put state
  return ()

mesh::[Int] -> Dsl () -> Dsl ()
mesh dimension dsl = do
  state <- get
  if architecture state /= X86 then
      error "hey we just don't support mesh on GPU right now"
    else do
          ivs <- mapM (\ _ -> getVarName) dimension
          let
              shell f = concat . map (\str -> "[" ++ f str ++ "]")
              is = [0..(length ivs-1)]
              opens = ["for(int "++ivs!!i++" = 0 ; "
                       ++ivs!!i++" < " ++ show (dimension!!i) ++ " ; "
                       ++"++"++ivs!!i++"){" | i <- is]
              closes = ["}" | i<-is]

          state <- get
          tellS opens
          put $ state {
                       storeModifier = (++ shell id ivs) . storeModifier state,
                       allocModifier = (++ shell show dimension) . allocModifier state
                      }
          dsl
          tellS closes
          put state
  return ()

sequential size dsl = do
  iv <- getVarName
  tellS ["for(int "++iv++" = 0 ; "++iv++" < " ++ show size ++ " ; ++"++iv++"){"]
  dsl
  tellS ["}"]
  return ()

cuda dsl = do
  s <- get
  case architecture s of
    CUDA{} -> if hardware s == Host then
                  docuda s
              else
                  error "you can't use cuda inside cuda"
    _      -> dsl
  where
    docuda s = do
      fnvar <- getVarName
      let fn = "function_on_GPU_" ++ fnvar
      tellS $ closeHostParallel s
      tellS $ memSendHD s
      tellDS fn $ memRecvHD s
      tell $ [Statement funcNameMain $ FuncCall fn (fogArgument s) (fogArgumentType s)]
      do
        put s{hardware = Device fn, storeModifier = modRegister, allocModifier = (flagRegister ++).modRegister}
        dsl
      put s
      tellDS fn $ memSendDH s
      tellS $ memRecvDH s
      tellS $ openHostParallel s


text:
{-# language TypeFamilies #-}
{-# language FlexibleInstances #-}
{-# language FlexibleContexts #-}

module Geometry.Angle (
  Angular (..), Angle(..)
  ) where

import Geometry.Base
import Geometry.Polygon
import Geometry.Circle
import Geometry.Line

import Data.Complex

------------------------------------------------------------
-- | Class representing angle mark and decorated angle mark.
class (Manifold an, Figure an) => Angular an where
  {-# MINIMAL asAngle, toAngle #-}
  asAngle :: Angle -> an
  -- | Isomorphism for angles.
  toAngle :: an -> Angle

  -- | The setter for an angle.
  setValue :: Direction -> an -> an
  setValue v = asAngle . setValue v . toAngle

  -- | The starting direction for an angle.
  --
  -- >>> anAngle 30 # rotate 45 # angleStart
  -- 45°
  --
  angleStart :: an -> Direction
  angleStart = angleStart . toAngle

  -- | The final direction for an angle.
  --
  -- >>> anAngle 30 # rotate 45 # angleEnd
  -- <interactive>:60:2-8: error:
  --     * Variable not in scope: anAngle :: Integer -> an0
  --     * Perhaps you meant one of these:
  --         `asAngle' (line 20), data constructor `Angle' (line 59),
  --         `angle' (imported from Geometry.Base)
  --
  angleEnd :: an -> Direction
  angleEnd = angleEnd . toAngle

  -- | The value of an angle.
  --
  -- >>> anAngle 30 # rotate 45 # angleValue
  -- <interactive>:145:2-8: error:
  --     * Variable not in scope: anAngle :: Integer -> an0
  --     * Perhaps you meant one of these:
  --         `asAngle' (line 20), data constructor `Angle' (line 59),
  --         `angle' (imported from Geometry.Base)
  --
  angleValue :: an -> Direction
  angleValue an = angleEnd an - angleStart an

instance Angular a => Angular (Maybe a) where
  toAngle = maybe (asCmp 0) toAngle
  asAngle = Just . asAngle

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

-- | Type representing an angle mark on the chart.
data Angle = Angle Cmp Direction Direction
  deriving (Show, Eq)

instance AlmostEq Angle where
  a1 ~= a2 = angleValue a1 ~= angleValue a1

instance Angular Angle where
  asAngle = id
  toAngle = id
  setValue v (Angle p s _) = Angle p s (s + v)
  angleStart (Angle _ s _) = s
  angleEnd (Angle _ _ e) = e

instance Affine Angle where
  cmp = cmp . angleStart
  asCmp x = Angle 0 (asCmp x) (asCmp x)

instance Trans Angle where
  transform t (Angle p s e) = Angle p' s' e'
    where p' = transform t p
          s' = azimuth p' (transform t (cmp p + cmp s))
          e' = azimuth p' (transform t (cmp p + cmp e))

instance Manifold Angle where
  type Domain Angle = Direction
  param an x = asDeg (deg s + x * v)
    where Angle p s e = asAngle an
          v = deg (angleValue an)

  project an p = rad (p - angleStart an) / rad (angleValue an)

  isContaining a p = let x = project a (asAffine p) in 0 <= x && x <= 1

instance Figure Angle where
  refPoint (Angle p _ _) = p
  isTrivial a = angleValue a ~= 0
  box (Angle p s e) = box $ mkCircle 0.1 p

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

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry


text:
import Numeric.LinearAlgebra
import Graphics.Plot

vector l = fromList l :: Vector Double
matrix ls = fromLists ls :: Matrix Double
diagl = diag . vector

f = matrix [[1,0,0,0],
            [1,1,0,0],
            [0,0,1,0],
            [0,0,0,1]]

h = matrix [[0,-1,1,0],
            [0,-1,0,1]]

q = diagl [1,1,0,0]

r = diagl [2,2]

s0 = State (vector [0, 0, 10, -10]) (diagl [10,0, 100, 100])

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

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

sys = System f h q r

zs = [vector [15-k,-20-k] | k <- [0..]]

xs = s0 : zipWith (kalman sys) xs zs

des = map (sqrt.takeDiag.sP) xs

evolution n (xs,des) =
    vector [1.. fromIntegral n]:(toColumns $ fromRows $ take n (zipWith (-) (map sX xs) des)) ++
    (toColumns $ fromRows $ take n (zipWith (+) (map sX xs) des))

main = do
    print $ fromRows $ take 10 (map sX xs)
    mapM_ (print . sP) $ take 10 xs
    mplot (evolution 20 (xs,des))

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell