text:
{-# LANGUAGE BangPatterns,
             ScopedTypeVariables,
             RecordWildCards,
             FlexibleContexts,
             TypeFamilies,
             DeriveGeneric #-}
module Main where

import Data.IDX ( decodeIDXFile, idxIntContent, IDXData )
{------------------------}
import Numeric.LinearAlgebra as NL ( Vector, fromList )
{------------------------}
import Data.Vector.Unboxed as DTU ( toList )

import Prelude hiding (readFile)
import Neuro
    ( Samples, Sample, createNetwork, tanh', saveNetwork, trainNTimes )

main :: IO ()
main = do
  n <- createNetwork 784 [64] 10
  samples <- importTrain
  {-----------------------------}
  let n' = trainNTimes 50 5.0 tanh tanh' n samples
  saveNetwork "smartNet5.nn" n'

importTrain :: IO (Samples Double)
importTrain = do
    Just idxTrain  <- decodeIDXFile "train-images.idx3-ubyte" -- image
    Just idxResult <- decodeIDXFile "train-labels.idx1-ubyte" -- result
    return $ samples idxTrain idxResult

samples :: IDXData -> IDXData -> [Sample Double]
samples idxTrain idxResult = Prelude.zip (image idxTrain) (result idxResult) :: [Sample Double]

image  = matrix2x2 . _0_1to0_255 . DTU.toList . idxIntContent
result = unitar . DTU.toList . idxIntContent

unitar :: [Int] -> [NL.Vector Double]
unitar []     = []
unitar (0:xs) = NL.fromList [1,-1,-1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (1:xs) = NL.fromList [-1,1,-1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (2:xs) = NL.fromList [-1,-1,1,-1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (3:xs) = NL.fromList [-1,-1,-1,1,-1,-1,-1,-1,-1,-1] : unitar xs
unitar (4:xs) = NL.fromList [-1,-1,-1,-1,1,-1,-1,-1,-1,-1] : unitar xs
unitar (5:xs) = NL.fromList [-1,-1,-1,-1,-1,1,-1,-1,-1,-1] : unitar xs
unitar (6:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,1,-1,-1,-1] : unitar xs
unitar (7:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,1,-1,-1] : unitar xs
unitar (8:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,-1,1,-1] : unitar xs
unitar (9:xs) = NL.fromList [-1,-1,-1,-1,-1,-1,-1,-1,-1,1] : unitar xs

_0_1to0_255 :: [Int] -> [Double]
_0_1to0_255 [] = []
_0_1to0_255 (x:xs) = (fromIntegral x / 255) : _0_1to0_255 xs

matrix2x2 :: [Double] -> [NL.Vector Double]
matrix2x2 [] = []
matrix2x2 xs = NL.fromList (Prelude.take 784 xs) : matrix2x2 (Prelude.drop 784 xs)

ext: hs


text:
{-# LANGUAGE PolyKinds #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE PackageImports #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE DefaultSignatures #-}
module Data.Roles
  ( Representational(rep)
  , Phantom(phantom)
  ) where

import Control.Applicative
import Control.Monad.ST
import Data.Complex
import Data.Monoid
import Data.Proxy
import Data.Type.Coercion
import Data.Coerce
import Unsafe.Coerce
import Data.Map
import Data.IntMap

class Representational (t :: k1 -> k2) where
  -- | An argument is representational if you can lift a coercion of the argument into one of the whole
  rep :: Coercion a b -> Coercion (t a) (t b)
  default rep :: Phantom t => Coercion a b -> Coercion (t a) (t b)
  rep _ = phantom

class Representational t => Phantom (t :: k1 -> k2) where
  -- | An argument is phantom if you can 'coerce' the whole ignoring the argument
  phantom :: Coercion (t a) (t b)
  default phantom :: Coercible (t a) (t b) => Coercion (t a) (t b)
  phantom = Coercion

-- * Data.Proxy

instance Representational Proxy
instance Phantom Proxy

-- * Const

instance Representational Const where rep Coercion = Coercion
instance Representational (Const a)
instance Phantom (Const a)

-- * Data.Type.Coercion

instance Representational Coercion     where rep = unsafeCoerce
instance Representational (Coercion a) where rep Coercion = Coercion

-- * Prelude

instance Representational (->)       where rep Coercion = Coercion
instance Representational ((->) a)   where rep Coercion = Coercion

instance Representational Either     where rep Coercion = Coercion
instance Representational (Either a) where rep Coercion = Coercion

instance Representational (,)     where rep Coercion = Coercion
instance Representational ((,) a) where rep Coercion = Coercion

instance Representational (,,)       where rep Coercion = Coercion
instance Representational ((,,) a)   where rep Coercion = Coercion
instance Representational ((,,) a b) where rep Coercion = Coercion

instance Representational (,,,)         where rep Coercion = Coercion
instance Representational ((,,,) a)     where rep Coercion = Coercion
instance Representational ((,,,) a b)   where rep Coercion = Coercion
instance Representational ((,,,) a b c) where rep Coercion = Coercion

instance Representational (,,,,)           where rep Coercion = Coercion
instance Representational ((,,,,) a)       where rep Coercion = Coercion
instance Representational ((,,,,) a b)     where rep Coercion = Coercion
instance Representational ((,,,,) a b c)   where rep Coercion = Coercion
instance Representational ((,,,,) a b c d) where rep Coercion = Coercion

instance Representational (,,,,,)             where rep Coercion = Coercion
instance Representational ((,,,,,) a)         where rep Coercion = Coercion
instance Representational ((,,,,,) a b)       where rep Coercion = Coercion
instance Representational ((,,,,,) a b c)     where rep Coercion = Coercion
instance Representational ((,,,,,) a b c d)   where rep Coercion = Coercion
instance Representational ((,,,,,) a b c d e) where rep Coercion = Coercion

instance Representational []      where rep Coercion = Coercion
instance Representational Maybe   where rep Coercion = Coercion
instance Representational IO      where rep Coercion = Coercion
instance Representational (ST s)  where rep Coercion = Coercion

-- * containers
instance Representational (Map k) where rep Coercion = Coercion
instance Representational IntMap  where rep Coercion = Coercion

-- * Data.Complex

instance Representational Complex where rep Coercion = Coercion

-- * Data.Monoid

instance Representational Sum     where rep Coercion = Coercion
instance Representational Product where rep Coercion = Coercion
instance Representational Dual    where rep Coercion = Coercion
instance Representational Endo    where rep Coercion = Coercion
instance Representational First   where rep Coercion = Coercion
instance Representational Last    where rep Coercion = Coercion

max_stars_repo_name: phadej/roles

max_issues_repo_name: phadej/roles

max_forks_repo_name: phadej/roles


text:
{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE KindSignatures        #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE StandaloneDeriving    #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE UndecidableInstances  #-}
{-|
Module      : Grenade.Core.Shape
Description : Dependently typed shapes of data which are passed between layers of a network
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

-}
module Grenade.Core.Shape (
    Shape (..)
  , S (..)
  , Sing (..)

  , randomOfShape
  , fromStorable
  ) where

import           Control.DeepSeq (NFData (..))
import           Control.Monad.Random ( MonadRandom, getRandom )

import           Data.Proxy
import           Data.Serialize
import           Data.Singletons
import           Data.Singletons.TypeLits
import           Data.Vector.Storable ( Vector )
import qualified Data.Vector.Storable as V

#if MIN_VERSION_base(4,11,0)
import           GHC.TypeLits hiding (natVal)
#else
import           GHC.TypeLits
#endif

import qualified Numeric.LinearAlgebra.Static as H
import           Numeric.LinearAlgebra.Static
import qualified Numeric.LinearAlgebra as NLA

-- | The current shapes we accept.
--   at the moment this is just one, two, and three dimensional
--   Vectors/Matricies.
--
--   These are only used with DataKinds, as Kind `Shape`, with Types 'D1, 'D2, 'D3.
data Shape
  = D1 Nat
  -- ^ One dimensional vector
  | D2 Nat Nat
  -- ^ Two dimensional matrix. Row, Column.
  | D3 Nat Nat Nat
  -- ^ Three dimensional matrix. Row, Column, Channels.

-- | Concrete data structures for a Shape.
--
--   All shapes are held in contiguous memory.
--   3D is held in a matrix (usually row oriented) which has height depth * rows.
data S (n :: Shape) where
  S1D :: ( KnownNat len )
      => R len
      -> S ('D1 len)

  S2D :: ( KnownNat rows, KnownNat columns )
      => L rows columns
      -> S ('D2 rows columns)

  S3D :: ( KnownNat rows
         , KnownNat columns
         , KnownNat depth
         , KnownNat (rows * depth))
      => L (rows * depth) columns
      -> S ('D3 rows columns depth)

deriving instance Show (S n)

-- Singleton instances.
--
-- These could probably be derived with template haskell, but this seems
-- clear and makes adding the KnownNat constraints simple.
-- We can also keep our code TH free, which is great.
data instance Sing (n :: Shape) where
  D1Sing :: Sing a -> Sing ('D1 a)
  D2Sing :: Sing a -> Sing b -> Sing ('D2 a b)
  D3Sing :: KnownNat (a * c) => Sing a -> Sing b -> Sing c -> Sing ('D3 a b c)

instance KnownNat a => SingI ('D1 a) where
  sing = D1Sing sing
instance (KnownNat a, KnownNat b) => SingI ('D2 a b) where
  sing = D2Sing sing sing
instance (KnownNat a, KnownNat b, KnownNat c, KnownNat (a * c)) => SingI ('D3 a b c) where
  sing = D3Sing sing sing sing

instance SingI x => Num (S x) where
  (+) = n2 (+)
  (-) = n2 (-)
  (*) = n2 (*)
  abs = n1 abs
  signum = n1 signum
  fromInteger x = nk (fromInteger x)

instance SingI x => Fractional (S x) where
  (/) = n2 (/)
  recip = n1 recip
  fromRational x = nk (fromRational x)

instance SingI x => Floating (S x) where
  pi = nk pi
  exp = n1 exp
  log = n1 log
  sqrt = n1 sqrt
  (**) = n2 (**)
  logBase = n2 logBase
  sin = n1 sin
  cos = n1 cos
  tan = n1 tan
  asin = n1 asin
  acos = n1 acos
  atan = n1 atan
  sinh = n1 sinh
  cosh = n1 cosh
  tanh = n1 tanh
  asinh = n1 asinh
  acosh = n1 acosh
  atanh = n1 atanh

--
-- I haven't made shapes strict, as sometimes they're not needed
-- (the last input gradient back for instance)
--
instance NFData (S x) where
  rnf (S1D x) = rnf x
  rnf (S2D x) = rnf x
  rnf (S3D x) = rnf x

-- | Generate random data of the desired shape
randomOfShape :: forall x m. ( MonadRandom m, SingI x ) => m (S x)
randomOfShape = do
  seed :: Int <- getRandom
  return $ case (sing :: Sing x) of
    D1Sing SNat ->
        S1D (randomVector  seed Uniform * 2 - 1)

    D2Sing SNat SNat ->
        S2D (uniformSample seed (-1) 1)

    D3Sing SNat SNat SNat ->
        S3D (uniformSample seed (-1) 1)

-- | Generate a shape from a Storable Vector.
--
--   Returns Nothing if the vector is of the wrong size.
fromStorable :: forall x. SingI x => Vector Double -> Maybe (S x)
fromStorable xs = case sing :: Sing x of
    D1Sing SNat ->
      S1D <$> H.create xs

    D2Sing SNat SNat ->
      S2D <$> mkL xs

    D3Sing SNat SNat SNat ->
      S3D <$> mkL xs
  where
    mkL :: forall rows columns. (KnownNat rows, KnownNat columns)
        => Vector Double -> Maybe (L rows columns)
    mkL v =
      let rows    = fromIntegral $ natVal (Proxy :: Proxy rows)
          columns = fromIntegral $ natVal (Proxy :: Proxy columns)
      in  if rows * columns == V.length v
             then H.create $ NLA.reshape columns v
             else Nothing

instance SingI x => Serialize (S x) where
  put i = (case i of
            (S1D x) -> putListOf put . NLA.toList . H.extract $ x
            (S2D x) -> putListOf put . NLA.toList . NLA.flatten . H.extract $ x
            (S3D x) -> putListOf put . NLA.toList . NLA.flatten . H.extract $ x
          ) :: PutM ()

  get = do
    Just i <- fromStorable . V.fromList <$> getListOf get
    return i

-- Helper function for creating the number instances
n1 :: ( forall a. Floating a => a -> a ) -> S x -> S x
n1 f (S1D x) = S1D (f x)
n1 f (S2D x) = S2D (f x)
n1 f (S3D x) = S3D (f x)

-- Helper function for creating the number instances
n2 :: ( forall a. Floating a => a -> a -> a ) -> S x -> S x -> S x
n2 f (S1D x) (S1D y) = S1D (f x y)
n2 f (S2D x) (S2D y) = S2D (f x y)
n2 f (S3D x) (S3D y) = S3D (f x y)

-- Helper function for creating the number instances
nk :: forall x. SingI x => Double -> S x
nk x = case (sing :: Sing x) of
  D1Sing SNat ->
    S1D (konst x)

  D2Sing SNat SNat ->
    S2D (konst x)

  D3Sing SNat SNat SNat ->
    S3D (konst x)

max_stars_repo_name: sebeaumont/grenade

max_issues_repo_name: sebeaumont/grenade

max_forks_repo_name: sebeaumont/grenade


text:
{-# LANGUAGE NoMonomorphismRestriction, ScopedTypeVariables, FlexibleContexts #-}

module Math.Probably.EMPCA where

import Numeric.LinearAlgebra hiding (orth)
import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.HMatrix ((#>), app, orth)
import Math.Probably.FoldingStats
import qualified Data.Vector.Storable as VS
import Foreign.Storable.Tuple ()
import Debug.Trace

import Data.Ord (comparing)
import Data.List (sortBy)

import Math.Probably.PCA
import Data.Random
import Unsafe.Coerce
import qualified Data.Binary as B
import Data.Word
import Control.Applicative
import Control.DeepSeq

data EmPcaBasis = EmPcaBasis {
   centering :: VS.Vector (Double,Double),
   emEvecs :: Mat
   } deriving (Show, Read)

instance B.Binary EmPcaBasis where
  put (EmPcaBasis c m) = do
      B.put (VS.length c)
      VS.mapM_ (\(c1,c2) -> B.put (toWord64 c1) >> B.put (toWord64 c2)) c
      B.put $ rows m
      B.put $ cols m
      VS.mapM_ (B.put . toWord64) $ flatten m
  get = do nc <- B.get
           c <- VS.replicateM nc $ (,) <$> fmap fromWord64 B.get <*> fmap fromWord64 B.get
           nrows <- B.get
           ncols <- B.get
           mflat <- VS.replicateM (nrows*ncols) $ fmap fromWord64 B.get
           return $ EmPcaBasis c $ reshape ncols mflat

instance NFData EmPcaBasis where
  rnf (EmPcaBasis x y) = rnf x `seq` rnf y

toWord64 :: a -> Word64
toWord64 = unsafeCoerce

fromWord64 :: Word64 -> a
fromWord64 = unsafeCoerce

emPca :: Int -> Int -> Maybe (Matrix Double) -> [Vec] -> RVar EmPcaBasis
emPca k iters mcinit vecs = do
  let meansds = findCentre vecs

      n = length vecs
      p = VS.length $ head vecs

      dat =  fromColumns $ map (centre meansds) vecs

      go c 0 = c
      go c iter =
        let x = inv (tr c <> c)<> tr c <> dat :: Matrix Double
            c1 = dat <> tr x <> inv (x<> tr x) :: Matrix Double
        in go (trace ("EM "++show iter ++ " C00="++show (c1!0!0)) c1) (iter-1)

  cinit <- case mcinit of
             Nothing -> fmap (p><k) $ sequence $ replicate (p*k) $ uniform 0 1
             Just c -> return c
  let cfinal = trdimit "cfinal" $ go cinit iters
      cOrth :: Matrix Double
      cOrth =  trdimit "cOrth" $ orth cfinal
      xevec = truePca $ trdimit "pcainput" $ tr cOrth <> dat
      evec = cOrth <>  xevec
--      truepca m = stat m

  --    (xevec, eval) = truepca (tr cOrth <> dat)
  return $ EmPcaBasis meansds evec

showC (m, iter) = "iter "++show iter++" C ="++dispf 2 m

applyEmPca :: EmPcaBasis -> Vec -> Vec
applyEmPca (EmPcaBasis cent vp) v =  tr vp <> centre cent v

findCentre :: [Vec] -> VS.Vector (Double,Double)
findCentre  = uncurry (VS.zipWith (,)) . runStat meanSDF

centre :: VS.Vector (Double,Double) -> Vec -> Vec
centre meansds v = VS.zipWith (\x (mn,_) -> x-mn) v meansds

trdims s m = trace (s++": "++show (rows m, cols m))

trdimit s m = trace (s++": "++show (rows m, cols m)) m

trdispit s m = trace (s++": "++dispf 3 m) m

trdisp s m = trace (s++": "++dispf 3 m)

truePca dataset =
  let --mn = trdispit "mn" $ asColumn (mean $ tr dataset_mn)
      --dataset = trdimit "dataset" $dataset_mn - mn
      cc = trdimit "cc" $ covN (tr dataset)
      (cdd,cvv) = eigSH cc
      ii = reverse $ map fst $ sortBy (comparing snd) $ zip [0..] $ VS.toList cdd
      evects = cvv ? ii
  in evects

max_stars_repo_name: glutamate/probably-baysig

max_issues_repo_name: silky/probably-baysig

max_forks_repo_name: silky/probably-baysig


text:
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}

module Lib
    ( TD(..),
      start,start',
      advances,advances',
      retreats,retreats',
      stepRandom,stepRandom',
      rhombus,
      diaTD
    ) where

import qualified Data.Map.Strict as M
import Data.Foldable
import Control.Monad
import Control.Monad.Random.Class
import Diagrams.Prelude hiding (conjugate,start)
import Data.Complex
import Data.Foldable

i = 0:+1
omega = (-1/2) :+ (sqrt 3/2)
omega' = conjugate omega

vec (x :+ y) = V2 x y
pt (x :+ y) = p2 (x,y)

data TD = TD {
  volume :: !Int,
  hXYtoZ :: M.Map (Int,Int) Int,
  hYZtoX :: M.Map (Int,Int) Int,
  hZXtoY :: M.Map (Int,Int) Int} deriving (Show,Eq)

start = TD 0 M.empty M.empty M.empty

startR :: Int -> Int -> M.Map (Int,Int) Int
startR m n = M.fromList [((i,j),0)|i<-[0..m-1],j<-[0..n-1]]

start' l m n = TD 0 (startR l m) (startR m n) (startR n l)

findC :: Int -> Int -> M.Map (Int,Int) Int -> Int
findC u v = M.findWithDefault 0 (u,v)

bound :: M.Map (Int,Int) Int -> Int
bound x = case (M.toDescList x) of
  [] -> 0
  ((a,_),_):_ -> a+1

adv :: Maybe Int -> Maybe Int
adv Nothing = Just 1
adv (Just x) = Just (x+1)

advance :: TD -> Int -> Int -> Int -> TD
advance (TD v xyZ yzX zxY) x y z =
  TD (v+1) (M.alter adv (x,y) xyZ) (M.alter adv (y,z) yzX) (M.alter adv (z,x) zxY)

advances :: TD -> [TD]
advances d@(TD v xyZ yzX zxY) = do
  x <- [0..bound xyZ]
  y <- [0..bound yzX]
  let z = findC x y xyZ
  guard $ x == findC y z yzX
  guard $ y == findC z x zxY
  return $ advance d x y z

advances' :: TD -> [TD]
advances' d@(TD v xyZ yzX zxY) = do
  ((x,y),z) <- M.toList xyZ
  guard $ Just x == M.lookup (y,z) yzX
  guard $ Just y == M.lookup (z,x) zxY
  return $ advance d x y z

rtr :: Maybe Int -> Maybe Int
rtr Nothing = Nothing -- shouldn't happen though
rtr (Just 1) = Nothing
rtr (Just x) = Just (x-1)

retreat :: TD -> Int -> Int -> Int -> TD
retreat (TD v xyZ yzX zxY) x y z =
  TD (v-1) (M.alter rtr (x,y) xyZ) (M.alter rtr (y,z) yzX) (M.alter rtr (z,x) zxY)

retreats :: TD -> [TD]
retreats d@(TD v xyZ yzX zxY) = do
  x <- [0..bound xyZ]
  y <- [0..bound yzX]
  let z = (findC x y xyZ) - 1
  guard $ z>=0
  guard $ x + 1 == findC y z yzX
  guard $ y + 1 == findC z x zxY
  return $ retreat d x y z

retreat' :: TD -> Int -> Int -> Int -> TD
retreat' (TD v xyZ yzX zxY) x y z =
  TD (v-1) (M.adjust s (x,y) xyZ) (M.adjust s (y,z) yzX) (M.adjust s (z,x) zxY)
  where s = subtract 1

retreats' :: TD -> [TD]
retreats' d@(TD v xyZ yzX zxY) = do
  ((x,y),z') <- M.toList xyZ
  let z = z'-1
  guard $ Just (x+1) == M.lookup (y,z) yzX
  guard $ Just (y+1) == M.lookup (z,x) zxY
  return $ retreat' d x y z

stepRandom d = uniform $ (advances d ++ retreats d)
stepRandom' d = uniform $ (advances' d ++ retreats' d)

rhombus :: Trail' Loop V2 Double
rhombus = fromVertices . map pt $ [0, omega, omega+omega', omega',0]

sh :: Int -> Int -> Int -> Complex Double
sh x y z = fromIntegral x*omega+fromIntegral y*omega'+fromIntegral z

rhombi = foldMap w . M.toList
  where w ((x,y),z) = strokeLoop rhombus # translate (vec $ sh x y z)

diaTD (TD _ xyZ yzX zxY) = (fc lightblue (rhombi  xyZ) <>
                            (rotateBy (1/3) . fc darkblue) (rhombi yzX) <>
                            (rotateBy (2/3) . fc blue) (rhombi zxY)) # lw thin # rotateBy (1/4) # pad 1.1

someFunc :: IO ()
someFunc = putStrLn "someFunc"


text:
module MNIST.DataSet where

import MNIST.Prelude

import qualified Codec.Compression.GZip  as GZip
import qualified Data.ByteString.Lazy    as L
import qualified Data.Vector.Generic     as G
import qualified Numeric.LinearAlgebra   as HM
import qualified Numeric.LinearAlgebra.Static  as SM

rows :: Int
rows = 28

cols :: Int
cols = 28

nPixels :: Int
nPixels = cols * rows

nLabels :: Int
nLabels = 10

_loadMNIST
    :: FilePath
    -> FilePath
    -> IO [(Int, UVector Int)]
_loadMNIST dataPath labelPath =
  runMaybeT (do
    i <- MaybeT . fmap (decodeIDX       . GZip.decompress) . L.readFile $ dataPath
    l <- MaybeT . fmap (decodeIDXLabels . GZip.decompress) . L.readFile $ labelPath
    d <- MaybeT . pure $ labeledIntData l i
    return d)
  >>= \case
    Just mnist -> return mnist
    Nothing    -> throwString $
      "couldn't read gzipped MNIST at data file " <> dataPath
      <> " and label file " <> labelPath

loadMNISTTf :: FilePath -> FilePath -> IO [(Int, Vector Int)]
loadMNISTTf dp lp = _loadMNIST dp lp
  >>= pure . fmap (identity *** G.convert)

trainingDataTf :: IO [(Int, Vector Int)]
trainingDataTf =
  loadMNISTTf
    "data/train-images-idx3-ubyte.gz"
    "data/train-labels-idx1-ubyte.gz"

testDataTf :: IO [(Int, Vector Int)]
testDataTf =
  loadMNISTTf
    "data/t10k-images-idx3-ubyte.gz"
    "data/t10k-labels-idx1-ubyte.gz"

-- ========================================================================= --

loadMNISTBp
    :: FilePath
    -> FilePath
    -> IO [(R 784, R 9)]
loadMNISTBp dp lp = _loadMNIST dp lp
  >>= pure . fmap ((fromJust . mkImage *** fromJust . mkLabel) . swap)

  where
    mkImage :: UVector Int -> Maybe (R 784)
    mkImage u = SM.create . G.convert . G.map (\i -> fromIntegral i / 255) $ u

    mkLabel :: Int -> Maybe (R 9)
    mkLabel n = SM.create $ HM.build 9 (fromIntegral . fromEnum . (== n) . round)

trainingDataBp :: IO [(R 784, R 9)]
trainingDataBp =
  loadMNISTBp
    "data/train-images-idx3-ubyte.gz"
    "data/train-labels-idx1-ubyte.gz"

testDataBp :: IO [(R 784, R 9)]
testDataBp =
  loadMNISTBp
    "data/t10k-images-idx3-ubyte.gz"
    "data/t10k-labels-idx1-ubyte.gz"

max_stars_repo_name: stites/haskell-mnist-benchmarks

max_issues_repo_name: stites/haskell-mnist-benchmarks

max_forks_repo_name: stites/haskell-mnist-benchmarks


text:
{-# language BangPatterns        #-}
{-# language LambdaCase          #-}
{-# language MagicHash           #-}
{-# language ScopedTypeVariables #-}

{-# OPTIONS_GHC -Wall #-}

module Streaming.FFT
  ( -- * streaming fft
    streamFFT

    -- * Types
  , Transform(..)
  , Window(..)
  ) where

import GHC.Classes (modInt#)
import Streaming.FFT.Types
import qualified Data.Complex as Complex
import qualified Data.Primitive.Contiguous as Contiguous
import qualified Data.Primitive.Contiguous.FFT as Contiguous
import qualified Streaming.Prelude as S
import Control.Monad.Trans.Class (lift)

intToDouble :: Int -> Double
{-# inline intToDouble #-}
intToDouble = fromIntegral

cissy :: Double -> Double -> Complex Double
cissy k n = Complex.cis (2 * pi * k / n)

windowSize :: Window -> Int
windowSize = \case
  Window64 -> 64
  Window128 -> 128
  Window256 -> 256
  Window512 -> 512
  Window1024 -> 1024

-- | Only safe when the second argument is not 0
unsafeMod :: Int -> Int -> Int
unsafeMod (I# x#) (I# y#) = I# (modInt# x# y#)
{-# inline unsafeMod #-} -- this should happen anyway. trust but verify.

-- | Compute the FFT of a previously computed FFT given a new sample.
--   This operation is done in-place.
--
--   /O(n)/
subFFT :: forall m. PrimMonad m
  => Complex Double -- ^ newest signal, x_{k+N}
  -> Transform m
  -> m ()
subFFT x_k_plus_N (Transform f1) = do
  n <- Contiguous.sizeMutable f1
  x_k <- Contiguous.read f1 0
  let go :: Int -> m ()
      go !f = if f < n
        then do {
            atIx <- Contiguous.read f1 f
          ; let expTerm = cissy (intToDouble $ (n + 1) * f) (intToDouble n)
                res = expTerm * (atIx + x_k_plus_N - x_k)
          ; Contiguous.write f1 f res
          ; go (f + 1)
        }
        else pure ()
  go 0

loadInitial :: forall m b. (PrimMonad m)
  => MutablePrimArray (PrimState m) (Complex Double)
  -> Window -- ^ window size
  -> Int -- ^ index
  -> Int -- ^ have we finished consuming the signal
  -> Stream (Of (Complex Double)) m b -- ^ first part of stream
  -> m (Stream (Of (Complex Double)) m b) -- ^ stream minus original signal
loadInitial marr sig !ix !untilSig stream =
  if (untilSig >= windowSize sig)
    then pure stream
    else do {
      S.next stream >>= \case
        Left _ -> pure stream
        Right (x,rest) -> if ix == 0
          then loadInitial marr sig (ix + 1) untilSig stream
          else do {
              Contiguous.write marr (unsafeMod (ix - 1 + untilSig) (windowSize sig)) x
            ; loadInitial marr sig (ix + 1) (untilSig + 1) rest
          }
    }

thereafter :: forall m b c. (PrimMonad m)
  => (Transform m -> m c)
  -> Window
  -> Int -- ^ have we filled the window size
  -> Transform m -- ^ transform
  -> Stream (Of (Complex Double)) m b
  -> Stream (Of c) m b
thereafter extract win !untilWin trans st =
  if (untilWin >= windowSize win)
    then thereafter extract win 0 trans st
    else do {
        e <- lift $ S.next st
      ; case e of {
            Left r -> pure r
          ; Right (x,rest) -> do {
                lift $ subFFT x trans
              ; info <- lift $ extract trans
              ; S.yield info
              ; thereafter extract win (untilWin + 1) trans rest
            }
        }
    }

streamFFT :: forall m b c. (PrimMonad m)
  => (Transform m -> m c)
  -> Window
  -> Stream (Of (Complex Double)) m b
  -> Stream (Of c) m b
streamFFT extract win stream = do
  -- Allocate the one array
  marr <- lift $ Contiguous.new (windowSize win)

  -- Grab the first signal from the stream
  streamMinusFirst <- lift $ loadInitial marr win 0 0 stream

  -- Get our first transform
  initialT :: Transform m <- lift $ do { Contiguous.mfft marr; pure $ Transform marr }

  -- Extract information from that transform
  initialInfo <- lift $ extract initialT

  -- Yield that information to the new stream
  S.yield initialInfo

  -- Now go
  thereafter extract win 0 initialT streamMinusFirst

max_stars_repo_name: chessai/streaming-fft

max_issues_repo_name: chessai/streaming-fft

max_forks_repo_name: chessai/streaming-fft


text:
module STCR2Z1T0PointSet where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:lenStr:initialScaleStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:histFilePath:pinwheelFlagStr:numIterationStr:writeSourceFlagStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      initialScale = read initialScaleStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z1T0PointSet"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z1T0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          sigma
          tao
          initialScale
          theta0Freqs
          thetaFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length theta0Freqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z1T0 <-
    computeUnboxedP $
    computeR2Z1T0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      thetaFreqs
      theta0Freqs
  plan <- makeR2Z1T0Plan emptyPlan arrR2Z1T0
  let n = 30
      m = round $ (fromIntegral n) * (sqrt 2) / 2
      a = 10
      b = -10
      c = 10
      a' = round $ (fromIntegral a) * (sqrt 2) / 2
      b' = round $ (fromIntegral b) * (sqrt 2) / 2
      c' = round $ (fromIntegral c) * (sqrt 2) / 2
      r = 30
      numTheta = 7
      deltaTheta = (1 * pi) / numTheta
      xs =
        ((L.map
            (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
            [ (r * cos (k * deltaTheta) + 0, r * sin (k * deltaTheta) + 0)
            | k <- [0 .. numTheta]
            ]) -- L.++
                -- (L.map
                --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
                --    [ ((r + 30) * cos (k * deltaTheta), (r + 30) * sin (k * deltaTheta))
                --    | k <- [0 .. numTheta - 1]
                --    ])
               -- L.++
                -- (L.map
                --    (\(i, j) -> R2S1RPPoint (round i, round j, 0, 1))
                --    [ ( (r + 15 * 2) * cos (k * deltaTheta)
                --      , (r + 15 * 2) * sin (k * deltaTheta))
                --    | k <- [0 .. numTheta - 1]
                --    ])
               -- L.++ [R2S1RPPoint (3, 5, 0, 1)]
         )
               -- ([R2S1RPPoint (i, i, 0, 1) | i <- [a',a' + b' .. c']] L.++
               --  [R2S1RPPoint (i, -i, 0, 1) | i <- [-a',-(a' + b') .. -c']] L.++
               --  [R2S1RPPoint (i, i, 0, 1) | i <- [-a',-(a' + b') .. -c']] L.++
               --  [R2S1RPPoint (i, -i, 0, 1) | i <- [a',a' + b' .. c']] L.++
               --  [R2S1RPPoint (i, 0, 0, 1) | i <- [a,a + b .. c]] L.++
               --  [R2S1RPPoint (0, i, 0, 1) | i <- [-a,-(a + b) .. -c]] L.++
               --  [R2S1RPPoint (i, 0, 0, 1) | i <- [-a,-(a + b) .. -c]] L.++
               --  [R2S1RPPoint (0, i, 0, 1) | i <- [a,a + b .. c]] L.++
               --  [R2S1RPPoint (3, 5, 0, 1)])
                -- [ R2S1RPPoint (n, 0, 0, 1)
                -- , R2S1RPPoint (0, n, 0, 1)
                -- , R2S1RPPoint (-n, 0, 0, 1)
                -- , R2S1RPPoint (0, -n, 0, 1)
                -- , R2S1RPPoint (m, m, 0, 1)
                -- , R2S1RPPoint (-m, m, 0, 1)
                -- , R2S1RPPoint (m, -m, 0, 1)
                -- , R2S1RPPoint (-m, -m, 0, 1)
                -- ]
  let bias = computeBiasR2T0 numPoint numPoint theta0Freqs xs
      eigenVec =
        computeInitialEigenVectorR2T0
          numPoint
          numPoint
          theta0Freqs
          thetaFreqs
          xs
  powerMethod1
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    arrR2Z1T0
    numIteration
    writeSourceFlag
    ""
    0.5
    bias
    eigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE AllowAmbiguousTypes #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE GADTs #-}
module Lib where

import Numeric
import Prelude hiding (and,or,(>))
import Control.Monad.RWS
import System.Process
import System.Exit
import Data.List (intercalate)
import qualified Data.Map as M
import qualified Data.Vector.Unboxed as V
import Statistics.Sample.KernelDensity
import Text.Read (readMaybe)

--------------------------------------------
-- Basic code-gen infrastructure

type a ⊸ b = a -> b
data Expr t where
  Expr :: String -> Expr t -- constant expression
  Lam :: (Expr a -> Expr t) -> Expr (a -> t) -- lambda
  App :: Expr (a -> t) -> Expr a -> Expr t -- application
  BinOp :: (String -> String -> String) -> Expr a -> Expr b -> Expr c
  UnOp :: (String -> String) -> Expr a -> Expr b
  If :: Prop -> Expr a -> Expr a -> Expr a
  -- Let :: P (Expr a) -> (Expr a ⊸ Expr b) -> Expr b
  Thunk :: P (Expr t) -> Expr t

class Representable a where
  constant :: a -> Expr a

instance Representable Bool where
  constant True = Expr "true"
  constant False = Expr "false"

instance Representable Float where
  constant = Expr . show

instance Fractional (Expr Float) where
  fromRational x = constant (fromRat x)
  (/) = BinOp (infixOp "/")

instance Num (Expr Float) where
  abs = UnOp (\x -> "Math.abs(" ++ x ++ ")")
  (-) = BinOp (infixOp "-")
  (+) = BinOp (infixOp "+")
  fromInteger x = constant (fromInteger x)
  (*) = BinOp (infixOp "*")

instance Floating (Expr Float) where
  sqrt = UnOp (\x -> "Math.sqrt(" ++ x ++ ")")
  exp = UnOp (\x -> "Math.exp(" ++ x ++ ")")
  log = UnOp (\x -> "Math.log(" ++ x ++ ")")

render :: Expr a -> P String
render = \case
  UnOp op x -> do
    x' <- render x
    return (op x')
  BinOp op x y -> do
    x' <- render x
    y' <- render y
    return (op x' y')
  Expr x -> return x
  App (Lam f) x -> render (f x)
  Lam f -> do
    x <- newVar
    body <- render (f (Expr x))
    return ("function( " ++ x ++ ") { return " ++ body ++ "; }")
  Thunk f -> do
    (x,body) <- censor (\_ -> []) $ listen f
    x' <- render x
    return ("function() {\n" ++ unlines body ++ " return " ++ x' ++ ";\n }")
  App x y -> do
    x' <- render x
    y' <- render y
    return ("(" ++ x' ++ ")(" ++ y' ++ ")")
  If cond x y -> do
    c' <- render cond
    x' <- render x
    y' <- render y
    return ("(" ++ c'  ++ "?" ++ x' ++ ":" ++ y' ++ ")")
  -- Let e f -> do
  --   x <- newVar
  --   e' <- render =<< e -- works only if f uses its argument at most once.
  --   fx <- render (f (Expr x))
  --   return ("var " ++ x ++ " = "++ e' ++";\n" ++ fx)

infixl #
(#) :: Expr (a -> b) -> Expr a -> Expr b
(#) = App

parens :: [Char] -> [Char]
parens x = "("++ x ++")"
infixOp :: [Char] -> [Char] -> [Char] -> [Char]
infixOp op x y = parens x ++ op ++ parens y
binFunc :: [Char] -> [Char] -> [Char] -> [Char]
binFunc f x y = f ++ "("++ x ++ "," ++ y ++")"

unFunc ::  [Char] -> [Char] -> [Char]
unFunc f x = f ++ "("++ x ++ ")"

newtype P a = P (RWS () [String] State a) deriving (Functor,Applicative,Monad,MonadState State, MonadWriter [String])

data State = State {nextVar :: Int}

logtrace :: String -> Expr a -> P ()
logtrace msg x = do
  x' <- render x
  emit ("console.log(" ++ show "LOG: " ++ "+" ++ show msg ++ "+ \" \" +" ++ x' ++ ");")

compileModel :: String -> P (Expr a) -> String
compileModel mainFunc m = unlines (z ++ [x])
  where (x,_,z) = runRWS p () (State {nextVar = 0})
        (P p) = render (UnOp (unFunc mainFunc) (Thunk m))

and :: Prop -> Prop -> Prop
and = BinOp (infixOp "&&")

(∧) :: Prop -> Prop -> Prop
(∧) = and

iff :: Prop -> Prop -> Prop
iff = BinOp (infixOp "==")

or :: Prop -> Prop -> Prop
or = BinOp (infixOp "||")

not' :: Prop -> Prop
not' = UnOp (\x -> "(!(" ++ x ++ "))")

(-->) :: Prop -> Prop -> Prop
p --> q = not' p `or` q

-- | Allocate a new variable
newVar :: P String
newVar = do
  n <- gets nextVar
  modify $ \State{..} -> State {nextVar = n+1, ..}
  return ("v" ++ show n)

emit :: String -> P ()
emit x = tell [x]

-----------------------------------
-- Types

type Vec = [Expr Float]
data Distrib a

type Ind = Vec
type Mat = [Vec]
-- type Vec = Expr Vector
type Prop = Expr Bool
type Pred = Ind -> Prop
type Measure = Ind -> Expr Float
type Adj = Vec
type AP = Measure
type CN = Ind -> Prop
type VP = Ind -> Prop
type NP = VP -> Prop
type Quant = CN -> NP

----------------------------------------------------
-- Compositional semantics

observe :: Prop -> P ()
observe = hyp

squared :: Num a => a -> a
squared x = x*x

observeEqual :: Expr Float -> Expr Float -> P ()
observeEqual x y = do
  f <- render (negate (squared (x-y)))
  emit ("factor(" ++ f ++ ");")

hyp :: Prop -> P ()
hyp x = do
  x' <- render x
  emit ("hyp(" ++ x' ++ ");")

-- | Sample new individual
newInd :: P Ind
newInd = newIndSuch []

-- | Sample new individual which satisfies some additional predicates
newIndSuch :: [Pred] -> P Ind
newIndSuch hs = do
  x <- newVector
  forM_ hs $ \h -> hyp (h x)
  return x

numberOfDimensions :: Int
numberOfDimensions = 2

newVector :: P Vec
newVector = mapM (uncurry sampleGaussian) (replicate numberOfDimensions (0,1))

newNormedVector :: P Vec
newNormedVector = do
  xs <- mapM (uncurry sampleGaussian) (replicate numberOfDimensions (0,1))
  return ((/ norm xs) <$> xs)

cosineDistance :: Vec -> Vec -> Expr Float
cosineDistance x y = dotProd x y / (norm x * norm y)

norm2 :: Vec -> Expr Float
norm2 x = dotProd x x

norm :: Vec -> Expr Float
norm = sqrt . norm2

newMatrix = mapM (\_ -> newVector) (replicate numberOfDimensions ())
newMatrix :: P Mat

satisfyAll :: [Ind] -> [Pred] -> P ()
satisfyAll xs ps = forM_ xs $ \x -> forM ps $ \p -> hyp (p x)

dotProd :: Vec -> Vec -> Expr Float
dotProd x y = sum (zipWith (*) x y)

vecMatProd :: Vec -> Mat -> Vec
vecMatProd v = map (dotProd v)

type Scalar = (Vec,Expr Float)

newScalarA :: P Scalar
newScalarA = do
  bias <- sampleGaussian 0 1
  v <- newNormedVector
  return (v,bias)

sampleGaussian :: Expr Float -> Expr Float -> P (Expr Float)
sampleGaussian mu sigma = do
  v <- newVar
  let m = Expr v
  mu' <- render mu
  sigma' <- render sigma
  emit ("var " ++ v ++ " = gaussian("  ++ mu' ++ "," ++ sigma' ++ ");")
  return m

newClass :: P Mat
newClass = newMatrix

forClass :: Mat -> Adj -> AP
forClass cl adj x = dotProd (vecMatProd adj cl) x

-- -- Alternative for scalar adjectives. We take the measure to be greater than that of a "random" element of the class.
-- forClassScalar :: Mat -> Adj -> AP
-- forClassScalar cl a x = Let (newIndSuch [isClass cl]) (\y -> adjAP a x - adjAP a y)

isClass :: Mat -> Pred
isClass clas x = dotProd (head clas) x > 0

adjAP :: Adj -> AP
adjAP = dotProd

gaussian :: Expr Float -> Expr Float -> Expr Float
gaussian mean stdDev = BinOp (binFunc "gaussian") mean stdDev

vague :: Float -> Measure -> Measure
vague vagueness m x = m x + gaussian 0 (Expr (show vagueness))

newMeasure :: P Measure
newMeasure = do
  (v,bias) <- newScalarA
  return (\x -> bias + dotProd v x)

positive :: Expr Float -> Prop
positive x = greaterThan x 0

newPred :: P (Ind -> Prop)
newPred = do
  m <- newMeasure
  return (\x -> positive (m x))

sample :: Expr (Distrib Bool) -> Expr Bool
sample d = Expr "sample" # d

bernouilli :: Expr Float -> Expr (Distrib Bool)
bernouilli = UnOp (\x -> "Bernoulli({p:" ++ x ++ "})")

-- TODO: there is a choice here.

-- this is an "intuitionistic" probably (it does not exclude the strong implication)
probablyInt :: Expr Float -> Prop -> Prop
probablyInt v x = sample (bernouilli v) --> x

-- this is a "definite" probably (it excludes the strong implication)
probablyDef ::  Expr Float -> Prop -> Prop
probablyDef v x = If (sample (bernouilli v)) x (not' x)

--  "a --> b"
-- and "if a then b else (not b)"
-- are not the same!

expectedPortion :: P Prop -> Expr Float
expectedPortion = UnOp (unFunc "expectedPortion") . Thunk

genQProb :: Expr Float -> Quant
genQProb prob cn vp = expectedPortion p > prob
  where p = do x <- newInd
               observe (cn x)
               return (vp x)

many :: Quant
many = genQProb 0.6

most :: Quant
most = genQProb 0.7

few :: Quant
few cn vp = genQProb 0.8 cn (not' . vp)

some :: Quant
some = genQProb 0.1

every :: Quant
every = genQProb 0.99

is :: Measure -> Pred
is m x =  (m x) > 0

more :: Measure -> Ind -> Ind -> Prop
more m x y = (m x) > (m y)

greaterThan :: Expr Float -> Expr Float -> Prop
greaterThan = BinOp (infixOp ">")

(>) :: Expr Float -> Expr Float -> Prop
(>) = greaterThan

-- equal :: Expr Float -> Expr Float -> Prop
-- equal x y = 0.1 > (abs (x - y))

equal :: Expr Float -> Expr Float -> Prop
equal x y = 0.1 > (abs (x - y))

disjoint :: Pred -> Pred -> Pred
disjoint p q x = not' (p x) `or` not' (q x)

subsective :: Adj -> Mat -> Pred
subsective a cl x = isClass cl x `and` is (forClass cl a) x

-- An alternative semantics for subsective scalars.

anything :: Pred
anything _ = Expr "true"

plot :: String -> [Double] -> IO ()
plot prefix xs = do
  let (xs',ys') = kde 64 (V.fromList xs)
      fname = prefix ++ ".svg"
  let ls =
        ["set terminal svg size 350,262 fname 'Verdana' enhanced background rgb 'white'",
         "set output '" ++ fname ++ "'",
         "set key off", -- no legend
          "$data << EOD"] ++
        [show x ++ " " ++ show y | (x,y) <- zip (V.toList xs') (V.toList ys')] ++
        ["EOD", "plot '$data' with lines"] -- filledcurve
  putStrLn "Plotting results..."
  (code,output,errors) <- readProcessWithExitCode "gnuplot" ["-p"] (unlines ls)
  case code of
    ExitFailure _ -> do
      putStrLn "Gnuplot failed with input:"
      putStrLn (unlines ls)
      putStrLn "errors:"
      putStrLn output
      putStrLn errors
    ExitSuccess ->
      putStrLn ("Plot output to " ++ fname)
  return ()

class KnownTyp a where
  isContinuous :: Bool

instance KnownTyp Float where
  isContinuous = True

instance KnownTyp Bool where
  isContinuous = False

parseValues :: [String] -> IO [Double]
parseValues [] = return []
parseValues (v:vs) = case readMaybe v of
  Nothing -> putStrLn v >> parseValues vs
  Just x -> (x:) <$> parseValues vs

run :: forall a. KnownTyp a => P (Expr a) -> IO ()
run model = do
  putStrLn "Creating model..."
  rts <- readFile "../Frontend/RTS.wppl"
  let m = compileModel mainName model
      funname = "modelFunction"
      mainName = if isContinuous @a then "mainContinuous" else "mainDiscrete"
      fname = funname ++ ".wppl"
  writeFile fname (intercalate "\n\n" [rts,m])
  putStrLn "Running webppl..."
  (code,output,errors) <- readProcessWithExitCode "webppl" [fname] ""
  case code of
    ExitFailure _ -> do
      putStrLn "Webppl failed with errors:"
      putStrLn output
      putStrLn errors
    ExitSuccess -> do
      putStrLn "Success!"
      case isContinuous @a of
        True -> do
            values <- parseValues (lines output)
            plot funname values
        False -> putStrLn output

max_stars_repo_name: GU-CLASP/CompositionalBayesianSemantics

max_issues_repo_name: GU-CLASP/CompositionalBayesianSemantics

max_forks_repo_name: GU-CLASP/CompositionalBayesianSemantics


text:
{-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-}

{- | The 'Numerical.HBLAS.BLAS.Level1' module provides a fully general
yet type safe Level1 BLAS API.

When in doubt about the semantics of an operation,
consult your system's BLAS api documentation, or just read the documentation
for
<https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm the Intel MKL BLAS distribution>

A few basic notes about how to invoke BLAS routines.

Many BLAS operations take one or more arguments of type 'Transpose'.
'Tranpose' has the following different constructors, which tell BLAS
routines what transformation to implicitly apply to an input matrix @mat@ with dimension @n x m@.

*  'NoTranspose' leaves the matrix @mat@ as is.

* 'Transpose' treats the @mat@ as being implicitly transposed, with dimension
    @m x n@. Entry @mat(i,j)@ being treated as actually being the entry
    @mat(j,i)@. For Real matrices this is also the matrix adjoint operation.
    ie @Tranpose(mat)(i,j)=mat(j,i)@

*  'ConjNoTranspose' will implicitly conjugate @mat@, which is a no op for Real ('Float' or 'Double') matrices, but for
'Complex Float' and 'Complex Double' matrices, a given matrix entry @mat(i,j)==x':+'y@
will be treated as actually being  @conjugate(mat)(i,j)=y':+'x@.

* 'ConjTranpose' will implicitly transpose and conjugate the input matrix.
ConjugateTranpose acts as matrix adjoint for both real and complex matrices.

The *gemm operations  work as follows (using 'sgemm' as an example):

* @'sgemm trLeft trRight alpha beta left right result'@, where @trLeft@ and @trRight@
are values of type 'Transpose' that respectively act on the matrices @left@ and @right@.

* the generalized matrix computation thusly formed can be viewed as being
@result = alpha * trLeft(left) * trRight(right) + beta * result@

the *gemv operations are akin to the *gemm operations, but with @right@ and @result@
being vectors rather than matrices.

the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@.
The 'MatUpLo' argument determines if the matrix should be treated as upper or
lower triangular and 'MatDiag' determines if the triangular solver should treat
the diagonal of the matrix as being all 1's or not.  A general pattern of invocation
would be @'strsv' matuplo  tranposeMatA  matdiag  matrixA  xVector@.
A key detail to note is that the input vector is ALSO the result vector,
ie 'strsv' and friends updates the vector place.

-}

module Numerical.HBLAS.BLAS.Level1(
        sasum
        ,dasum
        ,scasum
        ,dzasum

        ,saxpy
        ,daxpy
        ,caxpy
        ,zaxpy

        ,scopy
        ,dcopy
        ,ccopy
        ,zcopy

        ,sdot
        ,ddot
        --,sdsdot
        ,dsdot

        ,cdotu
        ,cdotc
        ,zdotu
        ,zdotc

        ,snrm2
        ,dnrm2
        ,scnrm2
        ,dznrm2

        ,srot
        ,drot

        --,srotg
        --,drotg

        ,srotm
        ,drotm

        --,srotmg
        --,drotmg

        ,sscal
        ,dscal
        ,cscal
        ,zscal
        ,csscal
        ,zdscal

        ,sswap
        ,dswap
        ,cswap
        ,zswap

        ,isamax
        ,idamax
        ,icamax
        ,izamax

            ) where

import Numerical.HBLAS.UtilsFFI
import Numerical.HBLAS.BLAS.FFI.Level1
import Numerical.HBLAS.BLAS.Internal.Level1
import Control.Monad.Primitive
import Data.Complex

sasum :: PrimMonad m => AsumFun Float (PrimState m) m Float
sasum = asumAbstraction "sasum" cblas_sasum_safe cblas_sasum_unsafe

dasum :: PrimMonad m => AsumFun Double (PrimState m) m Double
dasum = asumAbstraction "dasum" cblas_dasum_safe cblas_dasum_unsafe

scasum :: PrimMonad m => AsumFun (Complex Float) (PrimState m) m Float
scasum = asumAbstraction "scasum" cblas_scasum_safe cblas_scasum_unsafe

dzasum :: PrimMonad m => AsumFun (Complex Double) (PrimState m) m Double
dzasum = asumAbstraction "dzasum" cblas_dzasum_safe cblas_dzasum_unsafe

saxpy :: PrimMonad m => AxpyFun Float (PrimState m) m
saxpy = axpyAbstraction "saxpy" cblas_saxpy_safe cblas_saxpy_unsafe (\x f -> f x)

daxpy :: PrimMonad m => AxpyFun Double (PrimState m) m
daxpy = axpyAbstraction "daxpy" cblas_daxpy_safe cblas_daxpy_unsafe (\x f -> f x)

caxpy :: PrimMonad m => AxpyFun (Complex Float) (PrimState m) m
caxpy = axpyAbstraction "caxpy" cblas_caxpy_safe cblas_caxpy_unsafe withRStorable_

zaxpy :: PrimMonad m => AxpyFun (Complex Double) (PrimState m) m
zaxpy = axpyAbstraction "zaxpy" cblas_zaxpy_safe cblas_zaxpy_unsafe withRStorable_

scopy :: PrimMonad m => CopyFun Float (PrimState m) m
scopy = copyAbstraction "scopy" cblas_scopy_safe cblas_scopy_unsafe

dcopy :: PrimMonad m => CopyFun Double (PrimState m) m
dcopy = copyAbstraction "dcopy" cblas_dcopy_safe cblas_dcopy_unsafe

ccopy :: PrimMonad m => CopyFun (Complex Float) (PrimState m) m
ccopy = copyAbstraction "ccopy" cblas_ccopy_safe cblas_ccopy_unsafe

zcopy :: PrimMonad m => CopyFun (Complex Double) (PrimState m) m
zcopy = copyAbstraction "zcopy" cblas_zcopy_safe cblas_zcopy_unsafe

sdot :: PrimMonad m => NoScalarDotFun Float (PrimState m) m Float
sdot = noScalarDotAbstraction "sdot" cblas_sdot_safe cblas_sdot_unsafe

ddot :: PrimMonad m => NoScalarDotFun Double (PrimState m) m Double
ddot = noScalarDotAbstraction "ddot" cblas_ddot_safe cblas_ddot_unsafe

--sdsdot :: PrimMonad m => ScalarDotFun Float (PrimState m) m Float
--sdsdot = scalarDotAbstraction "sdsdot" cblas_sdsdot_safe cblas_sdsdot_unsafe withRStorable withRStorable

dsdot :: PrimMonad m => NoScalarDotFun Float (PrimState m) m Double
dsdot = noScalarDotAbstraction "dsdot" cblas_dsdot_safe cblas_dsdot_unsafe

cdotu :: PrimMonad m => ComplexDotFun (Complex Float) (PrimState m) m
cdotu = complexDotAbstraction "cdotu" cblas_cdotu_safe cblas_cdotu_unsafe

cdotc :: PrimMonad m => ComplexDotFun (Complex Float) (PrimState m) m
cdotc = complexDotAbstraction "cdotc" cblas_cdotc_safe cblas_cdotc_unsafe

zdotu :: PrimMonad m => ComplexDotFun (Complex Double) (PrimState m) m
zdotu = complexDotAbstraction "zdotu" cblas_zdotu_safe cblas_zdotu_unsafe

zdotc :: PrimMonad m => ComplexDotFun (Complex Double) (PrimState m) m
zdotc = complexDotAbstraction "zdotc" cblas_zdotc_safe cblas_zdotc_unsafe

snrm2 :: PrimMonad m => Nrm2Fun Float (PrimState m) m Float
snrm2 = norm2Abstraction "snrm2" cblas_snrm2_safe cblas_snrm2_unsafe

dnrm2 :: PrimMonad m => Nrm2Fun Double (PrimState m) m Double
dnrm2 = norm2Abstraction "dnrm2" cblas_dnrm2_safe cblas_dnrm2_unsafe

scnrm2 :: PrimMonad m => Nrm2Fun (Complex Float) (PrimState m) m Float
scnrm2 = norm2Abstraction "scnrm2" cblas_scnrm2_safe cblas_scnrm2_unsafe

dznrm2 :: PrimMonad m => Nrm2Fun (Complex Double) (PrimState m) m Double
dznrm2 = norm2Abstraction "dznrm2" cblas_dznrm2_safe cblas_dznrm2_unsafe

srot :: PrimMonad m => RotFun Float (PrimState m) m
srot = rotAbstraction "srot" cblas_srot_safe cblas_srot_unsafe

drot :: PrimMonad m => RotFun Double (PrimState m) m
drot = rotAbstraction "drot" cblas_drot_safe cblas_drot_unsafe

--srotg :: PrimMonad m => RotgFun Float (PrimState m) m
--srotg = rotgAbstraction "srotg" cblas_srotg_safe cblas_srotg_unsafe

--drotg :: PrimMonad m => RotgFun Double (PrimState m) m
--drotg = rotgAbstraction "drotg" cblas_drotg_safe cblas_drotg_unsafe

srotm :: PrimMonad m => RotmFun Float (PrimState m) m
srotm = rotmAbstraction "srotm" cblas_srotm_safe cblas_srotm_unsafe

drotm :: PrimMonad m => RotmFun Double (PrimState m) m
drotm = rotmAbstraction "drotm" cblas_drotm_safe cblas_drotm_unsafe

--srotmg :: PrimMonad m => RotmgFun Float (PrimState m) m
--srotmg = rotmgAbstraction "srotmg" cblas_srotmg_safe cblas_srotmg_unsafe

--drotmg :: PrimMonad m => RotmgFun Double (PrimState m) m
--drotmg = rotmgAbstraction "drotmg" cblas_drotmg_safe cblas_drotmg_unsafe

sscal :: PrimMonad m => ScalFun Float Float (PrimState m) m
sscal = scalAbstraction "sscal" cblas_sscal_safe cblas_sscal_unsafe (\x f -> f x )

dscal :: PrimMonad m => ScalFun Double Double (PrimState m) m
dscal = scalAbstraction "dscal" cblas_dscal_safe cblas_dscal_unsafe (\x f -> f x )

cscal :: PrimMonad m => ScalFun (Complex Float) (Complex Float) (PrimState m) m
cscal = scalAbstraction "cscal" cblas_cscal_safe cblas_cscal_unsafe withRStorable_

zscal :: PrimMonad m => ScalFun (Complex Double) (Complex Double) (PrimState m) m
zscal = scalAbstraction "zscal" cblas_zscal_safe cblas_zscal_unsafe withRStorable_

csscal :: PrimMonad m => ScalFun Float (Complex Float) (PrimState m) m
csscal = scalAbstraction "csscal" cblas_csscal_safe cblas_csscal_unsafe (\x f -> f x )

zdscal :: PrimMonad m => ScalFun Double (Complex Double) (PrimState m) m
zdscal = scalAbstraction "zdscal" cblas_zdscal_safe cblas_zdscal_unsafe (\x f -> f x )

sswap :: PrimMonad m => SwapFun Float (PrimState m) m
sswap = swapAbstraction "sswap" cblas_sswap_safe cblas_sswap_unsafe

dswap :: PrimMonad m => SwapFun Double (PrimState m) m
dswap = swapAbstraction "dswap" cblas_dswap_safe cblas_dswap_unsafe

cswap :: PrimMonad m => SwapFun (Complex Float) (PrimState m) m
cswap = swapAbstraction "cswap" cblas_cswap_safe cblas_cswap_unsafe

zswap :: PrimMonad m => SwapFun (Complex Double) (PrimState m) m
zswap = swapAbstraction "zswap" cblas_zswap_safe cblas_zswap_unsafe

isamax :: PrimMonad m => IamaxFun Float (PrimState m) m
isamax = iamaxAbstraction "isamax" cblas_isamax_safe cblas_isamax_unsafe

idamax :: PrimMonad m => IamaxFun Double (PrimState m) m
idamax = iamaxAbstraction "idamax" cblas_idamax_safe cblas_idamax_unsafe

icamax :: PrimMonad m => IamaxFun (Complex Float) (PrimState m) m
icamax = iamaxAbstraction "icamax" cblas_icamax_safe cblas_icamax_unsafe

izamax :: PrimMonad m => IamaxFun (Complex Double)(PrimState m) m
izamax = iamaxAbstraction "izamax" cblas_izamax_safe cblas_izamax_unsafe

{-
isamin :: PrimMonad m => IaminFun Float (PrimState m) m
isamin = iaminAbstraction "isamin" cblas_isamin_safe cblas_isamin_unsafe

idamin :: PrimMonad m => IaminFun Double (PrimState m) m
idamin = iaminAbstraction "idamin" cblas_idamin_safe cblas_idamin_unsafe

icamin :: PrimMonad m => IaminFun (Complex Float) (PrimState m) m
icamin = iaminAbstraction "icamin" cblas_icamin_safe cblas_icamin_unsafe

izamin :: PrimMonad m => IaminFun (Complex Double)(PrimState m) m
izamin = iaminAbstraction "izamin" cblas_izamin_safe cblas_izamin_unsafe
-}

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}

module HVX.Primitives
-- TODO(mh): Currently primitives that should generate implicit constraints are
-- unsupported. The should be supported shortily. (2014-06-04)
  ( apply
  , hadd
  , (+~)
  , hmul
  , (*~)
  , habs
  , neg
  , hlog
  , hexp
  , logsumexp
  , hmax
  , hmin
  , norm
  , berhu
  , huber
  , quadform
  , powBaseP0
  , powBaseP01
  , powBaseP1
  , powBaseP1InfEven
  , powBaseP1InfNotInt
  ) where

import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.HMatrix (eigSH')

import HVX.Internal.DCP
import HVX.Internal.Matrix
import HVX.Internal.Primitives
import HVX.Internal.Util

apply :: (ApplyVex vf mf ve me ~ vr, ValidVex vr)
  => Fun vf mf -> Expr ve me -> Expr vr (ApplyMon mf me)
apply = EFun

hadd :: (AddVex v1 v2 ~ v3, ValidVex v3)
  => Expr v1 m1 -> Expr v2 m2 -> Expr v3 (AddMon m1 m2)
hadd = EAdd

infixl 6 +~
(+~) :: (AddVex v1 v2 ~ v3, ValidVex v3)
  => Expr v1 m1 -> Expr v2 m2 -> Expr v3 (AddMon m1 m2)
(+~) = hadd

-- Constructors that enforce DCP constraints.
hmul :: (ApplyVex 'Affine 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Expr 'Affine 'Const -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
hmul (EConst a) e = apply (Mul a) e
hmul _ _ = error "the left argument of a multiply must be a constant"

infixl 7 *~
(*~) :: (ApplyVex 'Affine 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Expr 'Affine 'Const -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
(*~) (EConst a) e = apply (Mul a) e
(*~) _ _ = error "the left argument of a multiply must be a constant"

habs :: (ApplyVex 'Convex 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
habs = apply Abs

neg :: (ApplyVex 'Affine 'Noninc v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Noninc m1)
neg = apply Neg

hlog :: (ApplyVex 'Concave 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
hlog = apply Log

hexp :: (ApplyVex 'Convex 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
hexp e = apply Exp e

logsumexp :: (ApplyVex 'Convex 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
logsumexp = apply LogSumExp

hmax :: (ApplyVex 'Convex 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
hmax = apply Max

hmin :: (ApplyVex 'Concave 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
hmin = apply Min

norm :: (ApplyVex 'Convex 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
norm p
  | p == infinity = error "Internal: Infinity norm should become max . abs."
  | 1 <= p = apply (Norm p)
  | otherwise = error "Internal: Norm only supports p >= 1."

berhu :: (ApplyVex 'Convex 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
berhu m
  | 0 < m = apply (Berhu m)
  | otherwise = error "Internal: Berhu only supports m >= 0."

huber :: (ApplyVex 'Convex 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
huber m
  | 0 < m = apply (Huber m)
  | otherwise = error "Internal: Huber only supports m >= 0."

quadform :: (ApplyVex 'Convex 'Nonmon 'Affine m1 ~ v2, ValidVex v2)
  => Expr 'Affine 'Const -> Expr 'Affine m1 -> Expr v2 (ApplyMon 'Nonmon m1)
quadform (EConst a) e
  | rows a == cols a
    && fpequalsMat a (tr a)
    && 0 <= (maxElement.fst.eigSH' $ a) = apply (Quadform a) e -- we have checked that the matrix being passed in is Hermitian, so it's safe to use eigSH'.
  | otherwise = error "Matrices in quadratic forms must be positive semidefinite."
quadform _ _ = error "The matrix sandwitched by the quadratic form must be constant."

powBaseP0 :: (ApplyVex 'Affine 'Const v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Const m1)
powBaseP0 p
  | p `fpequals` 0 = apply PowBaseP0
  | otherwise = error "Internal: PowBaseP0 only supports p == 0."

powBaseP01 :: (ApplyVex 'Concave 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
powBaseP01 p
  | 0 < p && p < 1 = apply (PowBaseP01 p)
  | otherwise = error "Internal: PowBaseP01 only supports 0 < p < 1."

powBaseP1 :: (ApplyVex 'Affine 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
powBaseP1 p
  | p `fpequals` 1 = apply PowBaseP1
  | otherwise = error "Internal: PowBaseP1 only supports p == 1."

powBaseP1InfEven :: (ApplyVex 'Convex 'Nonmon v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nonmon m1)
powBaseP1InfEven p
  | 1 < p && isInteger p && even intP = apply (PowBaseP1InfEven intP)
  | otherwise = error "Internal: PowBaseP1InfEven only supports even p > 1."
  where intP = round p :: Integer

powBaseP1InfNotInt :: (ApplyVex 'Convex 'Nondec v1 m1 ~ v2, ValidVex v2)
  => Double -> Expr v1 m1 -> Expr v2 (ApplyMon 'Nondec m1)
powBaseP1InfNotInt p
  | 1 < p && not (isInteger p) = apply (PowBaseP1InfNotInt p)
  | otherwise = error "Internal: PowBaseP1InfNotInt only supports non integral p > 1."

max_stars_repo_name: chrisnc/hvx

max_issues_repo_name: chrisnc/hvx

max_forks_repo_name: chrisnc/hvx


text:
module Main (
    main,
    StatisticsSuite (..),
    runSuite
) where

import Prelude
import Statistics (runFold)
import qualified TextStatistics as TS
import Control.Applicative
import Data.Text.Lazy as T
import qualified Data.Text.Lazy.Encoding as E
import qualified Data.ByteString.Lazy as B

data StatisticsSuite = StatisticsSuite {
    wordCount :: Int,
    lineCount :: Int,
    averageWordLength :: Maybe Float,
    mostCommonLetter :: Maybe Char
} deriving (Eq, Ord, Show)

statisticsSuite :: TS.TextStatistics StatisticsSuite
statisticsSuite = StatisticsSuite <$> TS.wordCount <*> TS.lineCount <*> TS.averageWordLength <*> TS.mostCommonLetter

runSuite :: T.Text -> StatisticsSuite
runSuite = runFold foldl' statisticsSuite

main :: IO ()
main = do
    content <- fmap E.decodeUtf8 B.getContents
    print $ runSuite $ content

max_stars_repo_name: maciasoft/wc

max_issues_repo_name: maciasoft/wc

max_forks_repo_name: maciasoft/wc


text:
module Overture.Data ( module Data.Void
                     , module Data.Maybe
                     , module Data.Either
                     , module Data.These
                     , module Data.List
                     , module Data.Ord
                     , module Data.Monoid
                     , module Data.Foldable
                     , module Data.Traversable
                     , module Data.Complex
                     , module Data.Ix
                     -- * Renamed functions
                     , unit, (++)
                     -- * Generic equivalents
                     , map, replicate, drop, take, splitAt, length
                     -- * Replacements for needlessly partial functions
                     , head, tail, last, init, cycle, atIndex
                     , read, readEither
                     -- * Assorted utility functions
                     , enumerate, takeUntil2
                     , (??), justIf
                     , between
                     , rightIf, isLeft, isRight, eitherToMaybe, fromEither
                     ) where

import Data.Either
import Data.These
import Data.Void
import Data.Maybe ( maybe, fromMaybe, isNothing, isJust
                  , catMaybes, mapMaybe, maybeToList
                  )
import Data.List ( unfoldr, intersperse, intercalate, sort
                 , sortBy, groupBy, partition, stripPrefix
                 )
import Data.Ord (comparing)
import Data.Monoid (Monoid(..))
import Data.Foldable
import Data.Traversable
import Data.Complex (Complex(..), realPart, imagPart)
import Data.Ix

import Control.Applicative
import Text.ParserCombinators.ReadPrec
import Text.ParserCombinators.ReadP as ReadP hiding (between)
import Text.Read hiding (read)
import qualified Prelude
import qualified Data.List as List
import Overture.ExportPrelude
import Prelude ()

-- better names

-- | Alias for 'mempty'
unit :: (Monoid a) => a
unit = mempty

infixr 5 ++

-- | Alias for 'mappend'
(++) :: (Monoid a) => a -> a -> a
(++) = mappend

-- | A better name for 'fmap'.
map :: (Functor f) => (a -> b) -> (f a -> f b)
map = fmap

replicate :: (Integral z) => z -> a -> [a]
replicate = List.genericReplicate

drop :: (Integral z) => z -> [a] -> [a]
drop = List.genericDrop

take :: (Integral z) => z -> [a] -> [a]
take = List.genericTake

length :: (Integral z) => [a] -> z
length = List.genericLength

splitAt :: (Integral z) => z -> [a] -> ([a], [a])
splitAt = List.genericSplitAt

head :: [a] -> Maybe a
head []    = Nothing
head (x:_) = Just x

tail :: [a] -> Maybe [a]
tail []     = Nothing
tail (_:xs) = Just xs

last :: [a] -> Maybe a
last []     = Nothing
last [x]    = Just x
last (_:xs) = last xs

init :: [a] -> Maybe [a]
init []     = Nothing
init [_]    = Just []
init (x:xs) = (x:) <$> init xs

cycle :: [a] -> Maybe [a]
cycle [] = Nothing
cycle xs = Just $ Prelude.cycle xs

atIndex :: [a] -> Integer -> Maybe a
atIndex xs n = head $ drop n xs

readEither :: (Read a) => String -> Either String a
readEither s =
  case [ x | (x,"") <- readPrec_to_S read' minPrec s ] of
    [x] -> Right x
    []  -> Left "no parse"
    _   -> Left "ambiguous parse"
 where
  read' = readPrec >>= \x -> lift ReadP.skipSpaces >> return x

read :: (Read a) => String -> Maybe a
read s = case readEither s of
    Left  _ -> Nothing
    Right x -> Just x

-- | A list from 'minBound' to 'maxBound'.
enumerate :: (Enum a, Bounded a) => [a]
enumerate = [minBound .. maxBound]

-- | Takes a prefix of a list until two successive elements satisfy the
--   given predicate, keeping the first.
takeUntil2 :: (a -> a -> Bool) -> [a] -> [a]
takeUntil2 f (x1:x2:xs) | f x1 x2   = [x1]
                        | otherwise = x1:takeUntil2 f (x2:xs)
takeUntil2 _ xs = xs

-- | A synonym for 'fromMaybe' with inverted arguments. Right-associative to
--   allow chaining, e.g. @x ?? y ?? z@.
(??) :: Maybe a -> a -> a
(??) = flip fromMaybe

infixr 3 ??

-- | Keeps or discards values based on a predicate. Useful for unfolds.
justIf :: (a -> Bool) -> a -> Maybe a
justIf p x | p x 	   = Just x
           | otherwise = Nothing

-- | Encodes the result of a predicate using 'Either'.
rightIf :: (a -> Bool) -> (a -> Either a a)
rightIf p x = if p x then Right x else Left x

between :: (Ord a) => a -> a -> a -> Bool
between a z = (&&) <$> (> a) <*> (< z)

isLeft :: Either a b -> Bool
isLeft (Left _) = True
isLeft _        = False

isRight :: Either a b -> Bool
isRight (Right _) = True
isRight _         = False

eitherToMaybe :: Either a b -> Maybe b
eitherToMaybe = either (\_ -> Nothing) Just

fromEither :: (a -> b) -> Either a b -> b
fromEither f (Left l)  = f l
fromEither _ (Right r) = r

-- | Case analysis for 'Bool'.
bool :: r -> r -> Bool -> r
bool f _ False = f
bool _ t True  = t

-- some missing instances

instance Foldable ((,) a) where
    foldr f z (_, x) = f x z

instance Traversable ((,) a) where
    sequenceA (x, y) = (,) x <$> y

max_stars_repo_name: isomorphism/Overture

max_issues_repo_name: isomorphism/Overture

max_forks_repo_name: isomorphism/Overture


text:
{-# LANGUAGE FlexibleContexts #-}
-----------------------------------------------------------------------------
{- |
Module      :  Numeric.LinearAlgebra.Util
Copyright   :  (c) Alberto Ruiz 2013
License     :  GPL

Maintainer  :  Alberto Ruiz (aruiz at um dot es)
Stability   :  provisional

-}
-----------------------------------------------------------------------------

module Numeric.LinearAlgebra.Util(

    -- * Convenience functions
    size, disp,
    zeros, ones,
    diagl,
    row,
    col,
    (&),(!), (¦), (#),
    (?),(¿),
    rand, randn,
    cross,
    norm,
    unitary,
    mt,
    pairwiseD2,
    rowOuters,
    null1,
    null1sym,
    -- * Convolution
    -- ** 1D
    corr, conv, corrMin,
    -- ** 2D
    corr2, conv2, separable,
    -- * Tools for the Kronecker product
    --
    -- | (see A. Fusiello, A matter of notation: Several uses of the Kronecker product in
    --  3d computer vision, Pattern Recognition Letters 28 (15) (2007) 2127-2132)

    --
    -- | @`vec` (a \<> x \<> b) == ('trans' b ` 'kronecker' ` a) \<> 'vec' x@
    vec,
    vech,
    dup,
    vtrans
) where

import Numeric.Container
import Numeric.LinearAlgebra.Algorithms hiding (i)
import Numeric.Matrix()
import Numeric.Vector()

import System.Random(randomIO)
import Numeric.LinearAlgebra.Util.Convolution

disp :: Int -> Matrix Double -> IO ()
-- ^ show a matrix with given number of digits after the decimal point
disp n = putStrLn . dispf n

-- | pseudorandom matrix with uniform elements between 0 and 1
randm :: RandDist
     -> Int -- ^ rows
     -> Int -- ^ columns
     -> IO (Matrix Double)
randm d r c = do
    seed <- randomIO
    return (reshape c $ randomVector seed d (r*c))

-- | pseudorandom matrix with uniform elements between 0 and 1
rand :: Int -> Int -> IO (Matrix Double)
rand = randm Uniform

-- | pseudorandom matrix with normal elements
randn :: Int -> Int -> IO (Matrix Double)
randn = randm Gaussian

-- | create a real diagonal matrix from a list
diagl :: [Double] -> Matrix Double
diagl = diag . fromList

-- | a real matrix of zeros
zeros :: Int -- ^ rows
      -> Int -- ^ columns
      -> Matrix Double
zeros r c = konst 0 (r,c)

-- | a real matrix of ones
ones :: Int -- ^ rows
     -> Int -- ^ columns
     -> Matrix Double
ones r c = konst 1 (r,c)

-- | concatenation of real vectors
infixl 3 &
(&) :: Vector Double -> Vector Double -> Vector Double
a & b = join [a,b]

-- | horizontal concatenation of real matrices
infixl 3 !
(!) :: Matrix Double -> Matrix Double -> Matrix Double
a ! b = fromBlocks [[a,b]]

-- | (00A6) horizontal concatenation of real matrices
infixl 3 ¦
(¦) :: Matrix Double -> Matrix Double -> Matrix Double
a ¦ b = fromBlocks [[a,b]]

-- | vertical concatenation of real matrices
(#) :: Matrix Double -> Matrix Double -> Matrix Double
infixl 2 #
a # b = fromBlocks [[a],[b]]

-- | create a single row real matrix from a list
row :: [Double] -> Matrix Double
row = asRow . fromList

-- | create a single column real matrix from a list
col :: [Double] -> Matrix Double
col = asColumn . fromList

-- | extract selected rows
infixl 9 ?
(?) :: Element t => Matrix t -> [Int] -> Matrix t
(?) = flip extractRows

-- | (00BF) extract selected columns
infixl 9 ¿
(¿) :: Element t => Matrix t -> [Int] -> Matrix t
m ¿ ks = trans . extractRows ks . trans $ m

cross :: Vector Double -> Vector Double -> Vector Double
-- ^ cross product (for three-element real vectors)
cross x y | dim x == 3 && dim y == 3 = fromList [z1,z2,z3]
          | otherwise = error $ "cross ("++show x++") ("++show y++")"
  where
    [x1,x2,x3] = toList x
    [y1,y2,y3] = toList y
    z1 = x2*y3-x3*y2
    z2 = x3*y1-x1*y3
    z3 = x1*y2-x2*y1

norm :: Vector Double -> Double
-- ^ 2-norm of real vector
norm = pnorm PNorm2

-- | Obtains a vector in the same direction with 2-norm=1
unitary :: Vector Double -> Vector Double
unitary v = v / scalar (norm v)

-- | (rows &&& cols)
size :: Matrix t -> (Int, Int)
size m = (rows m, cols m)

-- | trans . inv
mt :: Matrix Double -> Matrix Double
mt = trans . inv

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

-- | Matrix of pairwise squared distances of row vectors
-- (using the matrix product trick in blog.smola.org)
pairwiseD2 :: Matrix Double -> Matrix Double -> Matrix Double
pairwiseD2 x y | ok = x2 `outer` oy + ox `outer` y2 - 2* x <> trans y
               | otherwise = error $ "pairwiseD2 with different number of columns: "
                                   ++ show (size x) ++ ", " ++ show (size y)
  where
    ox = one (rows x)
    oy = one (rows y)
    oc = one (cols x)
    one k = constant 1 k
    x2 = x * x <> oc
    y2 = y * y <> oc
    ok = cols x == cols y

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

-- | outer products of rows
rowOuters :: Matrix Double -> Matrix Double -> Matrix Double
rowOuters a b = a' * b'
  where
    a' = kronecker a (ones 1 (cols b))
    b' = kronecker (ones 1 (cols a)) b

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

-- | solution of overconstrained homogeneous linear system
null1 :: Matrix Double -> Vector Double
null1 = last . toColumns . snd . rightSV

-- | solution of overconstrained homogeneous symmetric linear system
null1sym :: Matrix Double -> Vector Double
null1sym = last . toColumns . snd . eigSH'

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

vec :: Element t => Matrix t -> Vector t
-- ^ stacking of columns
vec = flatten . trans

vech :: Element t => Matrix t -> Vector t
-- ^ half-vectorization (of the lower triangular part)
vech m = join . zipWith f [0..] . toColumns $ m
  where
    f k v = subVector k (dim v - k) v

dup :: (Num t, Num (Vector t), Element t) => Int -> Matrix t
-- ^ duplication matrix (@'dup' k \<> 'vech' m == 'vec' m@, for symmetric m of 'dim' k)
dup k = trans $ fromRows $ map f es
  where
    rs = zip [0..] (toRows (ident (k^(2::Int))))
    es = [(i,j) | j <- [0..k-1], i <- [0..k-1], i>=j ]
    f (i,j) | i == j = g (k*j + i)
            | otherwise = g (k*j + i) + g (k*i + j)
    g j = v
      where
        Just v = lookup j rs

vtrans :: Element t => Int -> Matrix t -> Matrix t
-- ^ generalized \"vector\" transposition: @'vtrans' 1 == 'trans'@, and @'vtrans' ('rows' m) m == 'asColumn' ('vec' m)@
vtrans p m | r == 0 = fromBlocks . map (map asColumn . takesV (replicate q p)) . toColumns $ m
           | otherwise = error $ "vtrans " ++ show p ++ " of matrix with " ++ show (rows m) ++ " rows"
  where
    (q,r) = divMod (rows m) p

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
{-# LANGUAGE BangPatterns, OverloadedStrings, RecordWildCards, RelaxedPolyRec, ViewPatterns #-}

module Network.HTTP.LoadTest.Report
    (
      reportBasic
    , reportEvents
    , reportFull
    , writeReport
    -- * Other reports
    , csvEvents
    -- * Helper functions
    , buildTime
    ) where

import Control.Monad (forM_)
import Criterion.Analysis (SampleAnalysis(..), OutlierEffect(..),
                           OutlierVariance(..))
import Data.Data (Data)
import Data.Function (on)
import Data.Maybe (fromMaybe)
import Data.Monoid (mappend, mconcat, mempty)
import Data.Text (Text)
import Data.Text.Buildable (build)
import Data.Text.Format (prec, shortest)
import Data.Text.Lazy.Builder (Builder)
import Data.Vector (Vector)
import Network.HTTP.LoadTest.Types (Analysis(..), Basic(..), Event(..),
                                    Summary(..), summEnd)
import Paths_pronk (getDataFileName)
import Prelude hiding (print)
import Statistics.Function (sort)
import Statistics.Resampling.Bootstrap (Estimate(..))
import Statistics.Sample.KernelDensity (kde)
import System.IO (Handle)
import System.IO.Unsafe (unsafePerformIO)
import Text.Hastache (MuType(..))
import Text.Hastache.Context (mkGenericContext)
import qualified Criterion.Report as R
import qualified Data.ByteString.Lazy as L
import qualified Data.HashMap.Strict as H
import qualified Data.List as List
import qualified Data.MeldableHeap as Q
import qualified Data.Text.Format as T
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Data.Vector as V
import qualified Text.Hastache as H

reportBasic :: Handle -> Analysis Basic -> IO ()
reportBasic h Analysis{..} = do
  let print a b = T.hprint h a b
  print "latency:\n" ()
  print "    mean:    {}\n" [time (mean latency)]
  print "    std dev: {}\n" [time (stdDev latency)]
  print "    99%:     {}\n    99.9%:   {}\n" (time latency99, time latency999)
  print "\nthroughput:  {}\n" [rate throughput]

reportFull :: (IO () -> IO ()) -> Handle -> Analysis SampleAnalysis -> IO ()
reportFull whenLoud h Analysis{..} = do
  let print a b = T.hprint h a b
  print "latency:\n" ()
  print "    mean:    {}\n" [time (estPoint (anMean latency))]
  whenLoud $ do
    print "      lower: {}\n" [time (estLowerBound (anMean latency))]
    print "      upper: {}\n" [time (estUpperBound (anMean latency))]
  print "    std dev: {}\n" [time (estPoint (anStdDev latency))]
  whenLoud $ do
    print "      lower: {}\n" [time (estLowerBound (anStdDev latency))]
    print "      upper: {}\n" [time (estUpperBound (anStdDev latency))]
  effect h (anOutlierVar latency)
  print "    99%:     {}\n    99.9%:   {}\n" (time latency99, time latency999)
  print "\nthroughput:  {}\n" [rate throughput]

time :: Double -> Builder
time = buildTime 4

rate :: Double -> Builder
rate r = prec 4 r `mappend` " req/sec"

buildTime :: Int -> Double -> Builder
buildTime precision t
     | t < 1e-3  = prec precision (t * 1e6) `mappend` " usec"
     | t < 1     = prec precision (t * 1e3) `mappend` " msec"
     | otherwise = prec precision t `mappend` " sec"

effect :: Handle -> OutlierVariance -> IO ()
effect h OutlierVariance{..} =
    case ovEffect of
      Unaffected -> return ()
      _ -> T.hprint h "    estimates {} affected by outliers ({}%)\n"
           (howMuch, T.fixed 1 (ovFraction * 100))
    where howMuch = case ovEffect of
                      Unaffected -> "not" :: Text
                      Slight     -> "slightly"
                      Moderate   -> "moderately"
                      Severe     -> "severely"

reportEvents :: Handle -> Vector Summary -> IO ()
reportEvents h sumv = do
  let evtMap = G.foldl' go H.empty . G.map summEvent $ sumv
      go m e = H.insertWith (+) (classify e) (1::Int) m
      classify Timeout          = 0
      classify HttpResponse{..} = respCode
  T.hprint h "responses:\n" ()
  forM_ (List.sort . H.toList $ evtMap) $ \(e,n) -> do
    let nameOf 0 = "timeout "
        nameOf k = "HTTP " `mappend` build k
    T.hprint h "    {} {}\n" (nameOf e, T.left 7 ' ' n)
  T.hprint h "\n" ()

csvEvents :: Vector Summary -> Builder
csvEvents sums = "start,elapsed,event\n" `mappend` G.foldr go mempty sums
  where
    firstStart = summStart (G.head sums)
    go Summary{..} b = mconcat [
                         shortest $ summStart - firstStart
                       , ","
                       , shortest summElapsed
                       , ","
                       , classify summEvent
                       , "\n"
                       ] `mappend` b
    classify Timeout          = "timeout"
    classify HttpResponse{..} = build respCode

-- | The path to the template and other files used for generating
-- reports.
templateDir :: FilePath
templateDir = unsafePerformIO $ getDataFileName "templates"
{-# NOINLINE templateDir #-}

writeReport :: (Data a) => FilePath -> Handle -> Double -> Analysis a -> IO ()
writeReport template h elapsed a@Analysis{..} = do
  let context "include" = MuLambdaM $
                          R.includeFile [templateDir, R.templateDir]
      context "elapsed"   = MuVariable elapsed
      context "latKdeTimes" = R.vector "x" latKdeTimes
      context "latKdePDF" = R.vector "x" latKdePDF
      context "latKde"    = R.vector2 "time" "pdf" latKdeTimes latKdePDF
      context "latValues" = MuList . map mkGenericContext . G.toList $ lats
      context "thrTimes" = R.vector "x" thrTimes
      context "thrValues" = R.vector "x" thrValues
      context "concTimes" = R.vector "x" . U.fromList $ map fstS conc
      context "concValues" = R.vector "x" . U.fromList $ map sndS conc
      context n = mkGenericContext a n
      (latKdeTimes,latKdePDF) = kde 128 . G.convert . G.map summElapsed $ latValues
      lats = G.map (\s -> s { summStart = summStart s - t }) latValues
          where t = summStart . G.head $ latValues
      (thrTimes,thrValues) = graphThroughput (min (G.length latValues) 50) elapsed latValues
      conc = graphConcurrency lats
  tpl <- R.loadTemplate [".",templateDir] template
  bs <- H.hastacheStr H.defaultConfig tpl context
  L.hPutStr h bs

data T = T (U.Vector Double) {-# UNPACK #-} !Double

-- | Compute a graph of throughput, requests completed per time
-- interval.
graphThroughput :: Int          -- ^ Number of time slices.
                -> Double       -- ^ Amount of time elapsed.
                -> V.Vector Summary -> (U.Vector Double, U.Vector Double)
graphThroughput slices elapsed sumv =
    (G.generate slices $ \i -> fromIntegral i * timeSlice,
     G.unfoldrN slices go (T endv 0))
  where go (T v i) = Just (fromIntegral (G.length a), T b j)
           where (a,b) = G.span (<=t) v
                 t = start + (j * timeSlice)
                 j = i+1
        timeSlice = elapsed / fromIntegral slices
        start = summStart . G.head $ sumv
        endv = G.convert . sort . G.map summEnd $ sumv

data S = S {
      fstS :: {-# UNPACK #-} !Double
    , sndS :: {-# UNPACK #-} !Int
    }

-- | Compute a graph of concurrency.
graphConcurrency :: V.Vector Summary -> [S]
graphConcurrency = scanl1 f . filter ((/=0) . sndS) . map (foldl1 (flip f)) .
                   List.groupBy ((==) `on` fstS) . go Q.empty . G.toList
  where
    f (S _ i) (S t j) = S t (i+j)
    go q es@(Summary{..}:xs)
        | summStart < t = S summStart 1 : go insQ xs
        | otherwise     = S t (-1)      : go delQ es
      where (t,delQ) = fromMaybe (1e300,q) $ Q.extractMin q
            insQ     = Q.insert (summStart+summElapsed) q
    go q _ = drain q
      where drain (Q.extractMin -> Just (t,q')) = S t (-1) : drain q'
            drain _ = []

max_stars_repo_name: fhartwig/pronk

max_issues_repo_name: liqd/pronk

max_forks_repo_name: bos/pronk


text:
module STCR2Z2T0S0EndModal where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           Image.Transform         (normalizeValueRange)
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:cutoffRadiusEndPointStr:cutoffRadiusStr:reversalFactorStr:cStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      cutoffRadiusEndPoint = read cutoffRadiusEndPointStr :: Int
      cutoffRadius = read cutoffRadiusStr :: Int
      reversalFactor = read reversalFactorStr :: Double
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0EndModal.noindex"
      a = 20
      b = 7
      c = read cStr :: Int
      endPointFilePath =
        folderPath </>
        (printf
           "EndPoint_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%d_%d_%d_%f.dat"
           numPoint
           (round thetaFreq :: Int)
           (round scaleFreq :: Int)
           (round maxScale :: Int)
           (round tao :: Int)
           cutoffRadiusEndPoint
           thetaSigma
           scaleSigma
           a
           b
           c
           reversalFactor)
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (cutoff cutoffRadius radialArr)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  arrR2Z2T0S0EndPoint <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (cutoff cutoffRadiusEndPoint radialArr)
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  let a' = round $ (fromIntegral a) * (sqrt 2) / 2
      b' = round $ (fromIntegral b) * (sqrt 2) / 2
      c' = round $ (fromIntegral c) * (sqrt 2) / 2
      xs =
        ([R2S1RPPoint (i, i, 0, 0) | i <- [a',a' + b' .. c']] L.++
         [R2S1RPPoint (i, -i, 0, 0) | i <- [-a',-(a' + b') .. -c']] L.++
         [R2S1RPPoint (i, i, 0, 0) | i <- [-a',-(a' + b') .. -c']] L.++
         [R2S1RPPoint (i, -i, 0, 0) | i <- [a',a' + b' .. c']] L.++
         [R2S1RPPoint (i, 0, 0, 0) | i <- [a,a + b .. c]] L.++
         [R2S1RPPoint (0, i, 0, 0) | i <- [-a,-(a + b) .. -c]] L.++
         [R2S1RPPoint (i, 0, 0, 0) | i <- [-a,-(a + b) .. -c]] L.++
         [R2S1RPPoint (0, i, 0, 0) | i <- [a,a + b .. c]])
      bias = computeBiasR2T0S0' numPoint numPoint theta0Freqs scale0Freqs xs
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          theta0Freqs
          scale0Freqs
          thetaFreqs
          scaleFreqs
          xs
  powerMethodR2Z2T0S0EndModal
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    numScale
    scaleFreqs
    scale0Freqs
    arrR2Z2T0S0EndPoint
    arrR2Z2T0S0
    numIteration
    writeSourceFlag
    (printf
       "_%d_%d_%d_%d_%d_%d_%.2f_%.2f_%f"
       numPoint
       (round thetaFreq :: Int)
       (round scaleFreq :: Int)
       (round maxScale :: Int)
       (round tao :: Int)
       cutoffRadiusEndPoint
       thetaSigma
       scaleSigma
       reversalFactor)
    reversalFactor
    bias
    eigenVec


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

module QUtil where

import Data.Complex
import Data.List

import Numeric.LinearAlgebra as L

import LinAlg
import BraKet
import Operators

-- Utility functions --
takeKet :: Element c => Int -> Ket c -> Ket c
takeKet n (Ket coeffs) = Ket $ L.subVector 0 n coeffs

-- Expected value of an operator O computed as Tr(rho <> O)
expectedValue :: MatOp -> [Complex Double] -> [Ket (Complex Double)] -> Complex Double
expectedValue op ps kets = trace $ op Prelude.<> densityMatrix ps kets

-- Commutator of two operators
commutator :: MatOp -> MatOp -> MatOp
commutator a b = (a Prelude.<> b) <~> (b Prelude.<> a)

-- normalization function
normalize :: Hilbert v => v -> v
normalize v = ((1.0 :+ 0.0) / (norm v :+ 0.0)) <**> v

-- sized vacuum Fock state
vacuum :: Int -> Ket (Complex Double)
vacuum n = Ket . L.fromList $ 1.0 : replicate (n-1) 0.0

-- n-photon Fock basis state
fockN :: Int -> Int -> Ket (Complex Double)
fockN m n = QUtil.normalize $ act (foldl1' (Prelude.<>) $ replicate n (create m)) (vacuum m)

-- Coherent state
coherent :: Int -> Complex Double -> Ket (Complex Double)
coherent m alpha = act (displacement m alpha) $ vacuum m

max_stars_repo_name: Jvanrhijn/optiqs

max_issues_repo_name: Jvanrhijn/optiqs

max_forks_repo_name: Jvanrhijn/optiqs


text:
module Test.StrictCheck.Consume
  ( Input
  , Inputs
  , Consume(..)
  , constructor
  , normalize
  , consumeTrivial
  , consumePrimitive
  ) where

import Test.QuickCheck
import Generics.SOP
import Generics.SOP.NS

import Test.StrictCheck.Internal.Inputs

import Data.Monoid    as Monoid
import Data.Semigroup as Semigroup
import Control.Applicative
import Control.Monad.State
import Control.Monad.Cont
import Control.Monad.Except
import Control.Monad.Identity
import Control.Monad.RWS
import Control.Monad.Reader
import Control.Monad.Writer
import Control.Monad.ST
import GHC.Conc
import Data.Functor.Product as Functor
import Data.Functor.Sum     as Functor
import Data.Functor.Compose as Functor
import Data.Complex

import Data.Foldable as Fold
import Data.List.NonEmpty (NonEmpty(..))
import Data.Tree     as Tree
import Data.Set      as Set
import Data.Map      as Map
import Data.Sequence as Seq
import Data.IntMap   as IntMap
import Data.IntSet   as IntSet

-------------------------------------------------------
-- The user interface for creating Consume instances --
-------------------------------------------------------

-- | Generate a tree of all possible ways to destruct the input value.
class Consume a where
  consume :: a -> Input
  default consume :: GConsume a => a -> Input
  consume = gConsume

-- | Reassemble pieces of input into a larger Input.
constructor :: Int -> [Input] -> Input
constructor n !is =
  Input (Variant (variant n)) is

-- | Use the CoArbitrary instance for a type to consume it. This should only be
-- used for "flat" types, i.e. those which contain no interesting substructure.
consumePrimitive :: CoArbitrary a => a -> Input
consumePrimitive !a =
  Input (Variant (coarbitrary a)) []

-- | If a type has no observable properties or substructure which can be used
-- to drive the randomization of output, consumption should merely evaluate a
-- value to weak-head normal form.
consumeTrivial :: a -> Input
consumeTrivial !_ =
  Input mempty []

-- | Fully normalize something which can be consumed
normalize :: Consume a => a -> ()
normalize (consume -> input) = go input
  where
    go (Input _ is) = Fold.foldr seq () (fmap go is)

--------------------------------------------
-- Deriving Consume instances generically --
--------------------------------------------

type GConsume a = (Generic a, All2 Consume (Code a))

gConsume :: GConsume a => a -> Input
gConsume !(from -> sop) =
  constructor (index_SOP sop)
  . hcollapse
  . hcliftA (Proxy :: Proxy Consume) (K . consume . unI)
  $ sop

---------------
-- Instances --
---------------

instance Consume (IO a)    where consume = consumeTrivial
instance Consume (STM a)   where consume = consumeTrivial
instance Consume (ST s a)  where consume = consumeTrivial
instance Consume (a -> b)  where consume = consumeTrivial
instance Consume (Proxy p) where consume = consumeTrivial

instance Consume Char     where consume = consumePrimitive
instance Consume Word     where consume = consumePrimitive
instance Consume Int      where consume = consumePrimitive
instance Consume Double   where consume = consumePrimitive
instance Consume Float    where consume = consumePrimitive
instance Consume Rational where consume = consumePrimitive
instance Consume Integer  where consume = consumePrimitive
instance (CoArbitrary a, RealFloat a) => Consume (Complex a) where
  consume = consumePrimitive

deriving newtype instance Consume a => Consume (Identity a)
deriving newtype instance Consume (RWST r w s m a)
deriving newtype instance Consume (ReaderT r m a)
deriving newtype instance Consume (StateT s m a)
deriving newtype instance Consume (m (a, w)) => Consume (WriterT w m a)
deriving newtype instance Consume (ContT r m a)
deriving newtype instance Consume (m (Either e a)) => Consume (ExceptT e m a)
deriving newtype instance Consume a => Consume (ZipList a)
deriving newtype instance Consume a => Consume (Monoid.First a)
deriving newtype instance Consume a => Consume (Monoid.Last a)
deriving newtype instance Consume a => Consume (Semigroup.First a)
deriving newtype instance Consume a => Consume (Semigroup.Last a)
deriving newtype instance Consume a => Consume (Min a)
deriving newtype instance Consume a => Consume (Max a)
deriving newtype instance Consume a => Consume (Monoid.Sum a)
deriving newtype instance Consume a => Consume (Monoid.Product a)
deriving newtype instance Consume (f (g a)) => Consume (Functor.Compose f g a)
deriving newtype instance Consume a => Consume (Dual a)
deriving newtype instance Consume a => Consume (Const a b)

instance Consume ()
instance Consume Bool
instance Consume Ordering
instance Consume a => Consume (Maybe a)
instance (Consume a, Consume b) => Consume (Either a b)
instance Consume a => Consume [a]

deriving newtype instance Consume a => Consume (I a)
deriving newtype instance Consume a => Consume (K a b)
-- TODO: instances for the rest of the SOP newtypes

instance (Consume (f a), Consume (g a)) => Consume (Functor.Sum f g a) where
  consume (InL a) = constructor 0 [consume a]
  consume (InR b) = constructor 1 [consume b]

instance (Consume (f a), Consume (g a)) => Consume (Functor.Product f g a) where
  consume (Pair a b) = constructor 0 [consume a, consume b]

instance (Consume a, Consume b) => Consume (Semigroup.Arg a b) where
  consume (Arg a b) = constructor 0 [consume a, consume b]

instance Consume a => Consume (NonEmpty a) where
  consume (a :| as) = constructor 0 [consume a, consume as]

instance Consume a => Consume (Tree a) where
  consume (Node a as) = constructor 0 [consume a, consume as]

instance Consume v => Consume (Map k v) where
  consume = constructor 0 . fmap (consume . snd) . Map.toList

consumeContainer :: (Consume a, Foldable t) => t a -> Input
consumeContainer = constructor 0 . fmap consume . Fold.toList

instance Consume v => Consume (Seq v)    where consume = consumeContainer
instance Consume v => Consume (Set v)    where consume = consumeContainer
instance Consume v => Consume (IntMap v) where consume = consumeContainer
instance Consume v => Consume IntSet where
  consume = consumeContainer . IntSet.toList

-- TODO: instances for the rest of Containers

instance (Consume a, Consume b) => Consume (a, b)
instance (Consume a, Consume b, Consume c) => Consume (a, b, c)
instance (Consume a, Consume b, Consume c, Consume d) => Consume (a, b, c, d)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e
         ) => Consume
  (a, b, c, d, e)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         ) => Consume
  (a, b, c, d, e, f)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g
         ) => Consume
  (a, b, c, d, e, f, g)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h
         ) => Consume
  (a, b, c, d, e, f, g, h)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i
         ) => Consume
  (a, b, c, d, e, f, g, h, i)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
          ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
         , Consume y
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y)
instance ( Consume a, Consume b, Consume c, Consume d, Consume e, Consume f
         , Consume g, Consume h, Consume i, Consume j, Consume k, Consume l
         , Consume m, Consume n, Consume o, Consume p, Consume q, Consume r
         , Consume s, Consume t, Consume u, Consume v, Consume w, Consume x
         , Consume y, Consume z
         ) => Consume
  (a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z)

max_stars_repo_name: nomeata/StrictCheck

max_issues_repo_name: nomeata/StrictCheck

max_forks_repo_name: nomeata/StrictCheck