text:
module SimpleParser.SimpleParser where

import System.IO
import Text.Parsec hiding ( spaces )
import Text.Parsec.String
import Control.Monad
import Numeric
import Data.Complex -- for complex number representation in 2.7
import Data.Ratio -- for rational numbers
import Data.Array -- could use haskell vectors instead, harder implementation
import Data.IORef
import Control.Monad.Error --deprecated but following the book :/

spaces :: Parser ()
spaces = skipMany1 space

type IOThrowsError = ErrorT LispError IO

data LispVal = Atom String
               | List [LispVal]
               | DottedList [LispVal] LispVal
               | Number Integer
               | String String
               | Bool Bool
               | Character Char
               | Float Double
               | Complex (Complex Double)
               | Rational Rational
               | Vector (Array Int LispVal)
               | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
               | Func { params :: [String], vararg :: Maybe String
                      , body :: [LispVal], closure :: Env}
               | IOFunc ([LispVal] -> IOThrowsError LispVal)
               | Port Handle

data LispError = NumArgs Integer [LispVal]
               | TypeMismatch String LispVal
               | Parser ParseError
               | BadSpecialForm String LispVal
               | NotFunction String String
               | UnboundVar String String
               | Default String

showError :: LispError -> String
showError (NumArgs expected found) = "Expected: " ++ show expected ++
  " args but found: " ++ unwordsList found
showError (TypeMismatch expected found) = "Type mismatch; Got: " ++
  show found ++ " but Expected: " ++ expected
showError (Parser error) = "Parse error at: " ++ show error
showError (UnboundVar message var) = message ++ " : " ++ show var
showError (BadSpecialForm message val) = message ++ " : " ++ show val
showError (NotFunction message func) = message ++ " : " ++ show func

instance Show LispError where show = showError
type ThrowsError = Either LispError
type Env = IORef [(String, IORef LispVal)]

parseAtom :: Parser LispVal
parseAtom =
  do
    first <- letter <|> symbol -- <|> is parsec choice operator
    rest <- many (letter <|> digit <|> symbol)
    let atom = first:rest
    return $ case atom of
      "#t" -> Bool True
      "#f" -> Bool False
      _    -> Atom atom

--------------------------- Exercise 2.1 ----------------------------------------
parseNumber' :: Parser LispVal
parseNumber' =
  do
    s <- many1 digit
    return (Number . read $ s)

--------------------------- Exercise 2.2 ----------------------------------------
parseNumber'' :: Parser LispVal
parseNumber'' = many1 digit >>= return . Number. read

--------------------------- Exercise 2.3 ----------------------------------------
specials :: Parser Char
specials = do char '\\'
              x <- oneOf "\\\"nrt"
              return $ case x of
                '\\' -> x
                '"'  -> x
                'n'  -> '\n'
                'r'  -> '\r'
                't'  -> '\t'

parseString :: Parser LispVal
parseString = do char '"'
                 x <- many $ specials <|> noneOf "\"\\"
                 char '"'
                 return $ String x
--------------------------- Exercise 2.4 ----------------------------------------
--Scheme defines Octal as #o, decimal as #d, hex as #h
symbol :: Parser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

-- so we can no longer parse bools with a prefixed '#'
parseBool :: Parser LispVal
parseBool =
  do
    char '#' --match a #
    (char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False))

parseNumber :: Parser LispVal
parseNumber = parseDecimal
  <|> parseSchemeDecimal
  <|> parseHex
  <|> parseOct
  <|> parseBin

parseDecimal :: Parser LispVal
parseDecimal = many1 digit >>= return . Number . read

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

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

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

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

numToList :: Int -> [Int]
numToList = map (read . (:[])) . show --probably slow, could use div and mod

oct2dig x = fst $ readOct x !! 0
hex2dig x = fst $ readHex x !! 0
bin2dig  = bin2dig' 0
bin2dig' digint "" = digint
bin2dig' digint (x:xs) = let old = 2 * digint + (if x == '0' then 0 else 1) in
                         bin2dig' old xs
--------------------------- Exercise 2.5 ----------------------------------------
parseCharacter :: Parser LispVal
--parseChar
             --try to parse a whole match on newline or space
             --or try to match any character not followed by an alpha numeric
             --return thst as a list
             -- now try to match newline or space return the corresponding value
             --if not a match then return the first character match by previous
             -- do
parseCharacter =
  do
    try $ string "#\\"
    x <- try (string "newline" <|> string "space")
         <|> do {str <- anyChar; notFollowedBy alphaNum; return [str]}
    return . Character $ case x of
      "newline" -> '\n'
      "space" -> ' '
      otherwise -> (x !! 0) --use of unsafe head

--------------------------- Exercise 2.6 ----------------------------------------
-- I keep struggling to find the appropriate informatin to match on in the
-- R5RS, I thought I had to implement exact/nonexact in addition to this...
parseFloat :: Parser LispVal
parseFloat =
  do
    x <- many1 digit --many  digits before a '.'
    char '.' -- match on the ','
    y <- many1 digit --more digitts after
    return $ Float . fst . head . readFloat $ (x ++ y) -- readFloat is a zipper

--------------------------- Exercise 2.7 ----------------------------------------
toDouble :: LispVal -> Double
toDouble (Float f) = realToFrac f
toDouble (Number n) = fromInteger n

parseComplexNumbers :: Parser LispVal
parseComplexNumbers =
  do
    a <- parseFloat <|> parseNumber
    char '+' --complex nums have pattern a+bi
    b <- parseFloat <|> parseNumber
    char 'i'
    return (Complex (toDouble a :+ toDouble b))

parseRationalNumbers :: Parser LispVal
parseRationalNumbers =
  do
    a <- many1 digit-- a rational num is denoted by a / b /= 0
    char '/'
    b <- many1 digit
    return (Rational ((read a) % (read b))) --didn't know that (%) existed

------------------------Begin section on recursive parsers-----------------------
parseList :: Parser LispVal
--sepBy parseExpr spaces seperates a string by spaces on parseExpr criteria
parseList = liftM List $ sepBy parseExpr spaces

parseDottedList :: Parser LispVal
parseDottedList =
  do
    head <- endBy parseExpr spaces
    tail <- char '.' >> spaces >> parseExpr
    return $ DottedList head tail

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

parseExpr :: Parser LispVal
parseExpr = parseAtom
  <|> parseString
  <|> try parseBool
  <|> try parseComplexNumbers
  <|> try parseRationalNumbers
  <|> parseNumber --try's are required because these start with a '#'
  <|> parseCharacter
  <|> parseQuasiQuoted
  <|> parseUnQuote
  <|> parseQuoted
  <|> try parseVector'
  <|> do char '(' --first bracker
         x <- try parseList <|> parseDottedList
         char ')'
         return x
--------------------------- Exercise 4.1 ----------------------------------------
parseQuasiQuoted :: Parser LispVal
parseQuasiQuoted =
  do
    char '`'
    x <- parseExpr
    return $ List [Atom "quasiquote", x]

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

--------------------------- Exercise 4.2 ----------------------------------------
parseVector :: Parser LispVal
parseVector =
  do
    values <- sepBy parseExpr spaces
    return $ Vector (listArray (0, (length values) - 1) values)

parseVector' :: Parser LispVal
parseVector' =
  do string "#("
     vector <- parseVector
     char ')'
     return vector

--------------------------- Exercise 4.3 ----------------------------------------
     --skip

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Float num) = show num
showVal (Complex num) = show num
showVal (Rational num) = show num
showVal (Vector arr) = unlines . map showVal . elems $ arr
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (List values) = "(" ++ unwordsList values ++ ")"
showVal (DottedList h t) = "(" ++ unwordsList h ++ " . " ++ showVal t
  ++ ")"
showVal (Character c) = [c]
showVal (PrimitiveFunc _) = "primitive"
showVal (Func {params = args, vararg = varargs, body = body, closure = env}) =
  "(lambda (" ++ unwords (map show args) ++
  (case varargs of
      Nothing -> ""
      Just arg -> " . " ++ arg) ++ ") ...)"
showVal (Port _) = "<IO port>"
showVal (IOFunc _) = "<IO primitive>"

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

instance Show LispVal where show = showVal

max_stars_repo_name: doyougnu/myLittleScheme

max_issues_repo_name: doyougnu/myLittleScheme

max_forks_repo_name: doyougnu/myLittleScheme


text:
-- |
-- Module    : Statistics.Distribution.Poisson.Internal
-- Copyright : (c) 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Internal code for the Poisson distribution.

module Statistics.Distribution.Poisson.Internal
    (
      probability, poissonEntropy
    ) where

import Data.List (unfoldr)
import Numeric.MathFunctions.Constants (m_sqrt_2_pi, m_tiny, m_epsilon)
import Numeric.SpecFunctions (logGamma, stirlingError {-, choose, logFactorial -})
import Numeric.SpecFunctions.Extra (bd0)

-- | An unchecked, non-integer-valued version of Loader's saddle point
-- algorithm.
probability :: Double -> Double -> Double
probability 0      0     = 1
probability 0      1     = 0
probability lambda x
  | isInfinite lambda    = 0
  | x < 0                = 0
  | x <= lambda * m_tiny = exp (-lambda)
  | lambda < x * m_tiny  = exp (-lambda + x * log lambda - logGamma (x+1))
  | otherwise            = exp (-(stirlingError x) - bd0 x lambda) /
                           (m_sqrt_2_pi * sqrt x)

-- -- | Compute entropy using Theorem 1 from "Sharp Bounds on the Entropy
-- -- of the Poisson Law".  This function is unused because 'directEntorpy'
-- -- is just as accurate and is faster by about a factor of 4.
-- alyThm1 :: Double -> Double
-- alyThm1 lambda =
--   sum (takeWhile (\x -> abs x >= m_epsilon * lll) alySeries) + lll
--   where lll = lambda * (1 - log lambda)
--         alySeries =
--           [ alyc k * exp (fromIntegral k * log lambda - logFactorial k)
--           | k <- [2..] ]

-- alyc :: Int -> Double
-- alyc k =
--   sum [ parity j * choose (k-1) j * log (fromIntegral j+1) | j <- [0..k-1] ]
--   where parity j
--           | even (k-j) = -1
--           | otherwise  = 1

-- | Returns [x, x^2, x^3, x^4, ...]
powers :: Double -> [Double]
powers x = unfoldr (\y -> Just (y*x,y*x)) 1

-- | Returns an upper bound according to theorem 2 of "Sharp Bounds on
-- the Entropy of the Poisson Law"
alyThm2Upper :: Double -> [Double] -> Double
alyThm2Upper lambda coefficients =
  1.4189385332046727 + 0.5 * log lambda +
  zipCoefficients lambda coefficients

-- | Returns the average of the upper and lower bounds according to
-- theorem 2.
alyThm2 :: Double -> [Double] -> [Double] -> Double
alyThm2 lambda upper lower =
  alyThm2Upper lambda upper + 0.5 * (zipCoefficients lambda lower)

zipCoefficients :: Double -> [Double] -> Double
zipCoefficients lambda coefficients =
  (sum $ map (uncurry (*)) (zip (powers $ recip lambda) coefficients))

-- Mathematica code deriving the coefficients below:
--
-- poissonMoment[0, s_] := 1
-- poissonMoment[1, s_] := 0
-- poissonMoment[k_, s_] :=
--   Sum[s * Binomial[k - 1, j] * poissonMoment[j, s], {j, 0, k - 2}]
--
-- upperSeries[m_]  :=
--  Distribute[Integrate[
--    Sum[(-1)^(j - 1) *
--      poissonMoment[j, \[Lambda]] / (j * (j - 1)* \[Lambda]^j),
--     {j, 3, 2 m - 1}],
--    \[Lambda]]]
--
-- lowerSeries[m_] :=
--  Distribute[Integrate[
--    poissonMoment[
--      2 m + 2, \[Lambda]] / ((2 m +
--         1)*\[Lambda]^(2 m + 2)), \[Lambda]]]
--
-- upperBound[m_] := upperSeries[m] + (Log[2*Pi*\[Lambda]] + 1)/2
--
-- lowerBound[m_] := upperBound[m] + lowerSeries[m]

upperCoefficients4 :: [Double]
upperCoefficients4 = [1/12, 1/24, -103/180, -13/40, -1/210]

lowerCoefficients4 :: [Double]
lowerCoefficients4 = [0,0,0, -105/4, -210, -2275/18, -167/21, -1/72]

upperCoefficients6 :: [Double]
upperCoefficients6 = [1/12, 1/24, 19/360, 9/80, -38827/2520,
                      -74855/1008, -73061/2520, -827/720, -1/990]

lowerCoefficients6 :: [Double]
lowerCoefficients6 = [0,0,0,0,0, -3465/2, -45045, -466235/4, -531916/9,
                      -56287/10, -629/11, -1/156]

upperCoefficients8 :: [Double]
upperCoefficients8 = [1/12, 1/24, 19/360, 9/80, 863/2520, 1375/1008,
                      -3023561/2520, -15174047/720, -231835511/5940,
                      -18927611/1320, -58315591/60060, -23641/3640,
                      -1/2730]

lowerCoefficients8 :: [Double]
lowerCoefficients8 = [0,0,0,0,0,0,0, -2027025/8, -15315300, -105252147,
                      -178127950, -343908565/4, -10929270, -3721149/14,
                      -7709/15, -1/272]

upperCoefficients10 :: [Double]
upperCoefficients10 = [1/12, 1/24, 19/360, 9,80, 863/2520, 1375/1008,
                       33953/5040, 57281/1440, -2271071617/11880,
                       -1483674219/176, -31714406276557/720720,
                       -7531072742237/131040, -1405507544003/65520,
                       -21001919627/10080, -1365808297/36720,
                       -26059/544, -1/5814]

lowerCoefficients10 :: [Double]
lowerCoefficients10 = [0,0,0,0,0,0,0,0,0,-130945815/2, -7638505875,
                       -438256243425/4, -435477637540, -3552526473925/6,
                       -857611717105/3, -545654955967/12, -5794690528/3,
                       -578334559/42, -699043/133, -1/420]

upperCoefficients12 :: [Double]
upperCoefficients12 = [1/12, 1/24, 19/360, 863/2520, 1375/1008,
                       33953/5040, 57281/1440, 3250433/11880,
                       378351/176, -37521922090657/720720,
                       -612415657466657/131040, -3476857538815223/65520,
                       -243882174660761/1440, -34160796727900637/183600,
                       -39453820646687/544, -750984629069237/81396,
                       -2934056300989/9576, -20394527513/12540,
                       -3829559/9240, -1/10626]

lowerCoefficients12 :: [Double]
lowerCoefficients12 = [0,0,0,0,0,0,0,0,0,0,0,
                       -105411381075/4, -5270569053750, -272908057767345/2,
                       -1051953238104769, -24557168490009155/8,
                       -3683261873403112, -5461918738302026/3,
                       -347362037754732, -2205885452434521/100,
                       -12237195698286/35, -16926981721/22,
                       -6710881/155, -1/600]

-- | Compute entropy directly from its definition. This is just as accurate
-- as 'alyThm1' for lambda <= 1 and is faster, but is slow for large lambda,
-- and produces some underestimation due to accumulation of floating point
-- error.
directEntropy :: Double -> Double
directEntropy lambda =
  negate . sum $
  takeWhile (< negate m_epsilon * lambda) $
  dropWhile (not . (< negate m_epsilon * lambda)) $
  [ let x = probability lambda k in x * log x | k <- [0..]]

-- | Compute the entropy of a Poisson distribution using the best available
-- method.
poissonEntropy :: Double -> Double
poissonEntropy lambda
  | lambda == 0 = 0
  | lambda <= 10 = directEntropy lambda
  | lambda <= 12 = alyThm2 lambda upperCoefficients4 lowerCoefficients4
  | lambda <= 18 = alyThm2 lambda upperCoefficients6 lowerCoefficients6
  | lambda <= 24 = alyThm2 lambda upperCoefficients8 lowerCoefficients8
  | lambda <= 30 = alyThm2 lambda upperCoefficients10 lowerCoefficients10
  | otherwise = alyThm2 lambda upperCoefficients12 lowerCoefficients12

ext: hs

max_stars_repo_name: vaerksted/statistics

max_issues_repo_name: vaerksted/statistics

max_forks_repo_name: vaerksted/statistics


text:
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE OverloadedStrings         #-}
{-# LANGUAGE PolyKinds                 #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# OPTIONS_GHC -fwarn-incomplete-patterns #-}

module Math.HMatrixUtils where

import qualified Knit.Effect.Logger            as KL
import qualified Knit.Report                   as K
import qualified Data.Text                     as T

import qualified Numeric.LinearAlgebra         as LA
import           Numeric.LinearAlgebra.Data     ( Matrix
                                                , R
                                                , Vector
                                                )

textSize :: (LA.Container c e, Show (LA.IndexOf c)) => c e -> T.Text
textSize = T.pack . show . LA.size

checkEqualVectors
  :: K.Member (KL.Logger KL.LogEntry) effs
  => T.Text
  -> T.Text
  -> Vector R
  -> Vector R
  -> K.Sem effs ()
checkEqualVectors nA nB vA vB = if LA.size vA == LA.size vB
  then return ()
  else
    KL.logLE KL.Error
    $  "Unequal vector length. length("
    <> nA
    <> ")="
    <> textSize vA
    <> " and length("
    <> nB
    <> ")="
    <> textSize vB

checkMatrixVector
  :: K.Member (KL.Logger KL.LogEntry) effs
  => T.Text
  -> T.Text
  -> Matrix R
  -> Vector R
  -> K.Sem effs ()
checkMatrixVector nA nB mA vB = if snd (LA.size mA) == LA.size vB
  then return ()
  else
    KL.logLE KL.Error
    $  "Bad matrix * vector lengths. dim("
    <> nA
    <> ")="
    <> textSize mA
    <> " and length("
    <> nB
    <> ")="
    <> textSize vB

checkVectorMatrix
  :: K.Member (KL.Logger KL.LogEntry) effs
  => T.Text
  -> T.Text
  -> Vector R
  -> Matrix R
  -> K.Sem effs ()
checkVectorMatrix nA nB vA mB = if LA.size vA == fst (LA.size mB)
  then return ()
  else
    KL.logLE KL.Error
    $  "Bad vector * matrix lengths. length("
    <> nA
    <> ")="
    <> textSize vA
    <> " and dim("
    <> nB
    <> ")="
    <> textSize mB

max_stars_repo_name: teto/Frames-utils

max_issues_repo_name: teto/Frames-utils

max_forks_repo_name: teto/Frames-utils


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

module BNN where

import System.Random.MWC (createSystemRandom)

--import Data.Aeson (ToJSON(toJSON), Value)
--import Data.Aeson (encode)
import qualified Data.ByteString.Lazy.Char8 as BL

import Data.List (sort, partition)
import qualified Data.Histogram as DH
import qualified Data.Histogram.Fill as DF
import qualified Data.Histogram.Bin.Bin2D as Bin2D
import Data.Vector.Unboxed.Base (Unbox)
import qualified Data.Vector.Unboxed as Vec

import Control.Monad (liftM2, replicateM, forM, forM_)
import Control.Monad.IO.Class (liftIO)
import Control.Monad.Bayes.Class
import Control.Monad.Bayes.Sampler
import Control.Monad.Bayes.Traced
import Control.Monad.Bayes.Weighted
import Control.Monad.Bayes.Inference.SMC as SMC
import Control.Monad.Bayes.Inference.RMSMC as RMSMC
import Control.Monad.Bayes.Sequential
import Control.Monad.Bayes.Population
import Control.Monad.Bayes.Traced.Static (Traced)
import Numeric.LinearAlgebra (Matrix, Vector, toList, vector, matrix, dot, (#>), (!), cmap, scalar, size)

import Numeric.Log

-- Plotting

histo2D :: (Foldable f, Unbox val, Num val) =>
    (Double, Int, Double)
    -> (Double, Int, Double)
    -> f (Double, Double)
    -> DH.Histogram (Bin2D.Bin2D DH.BinD DH.BinD) val
histo2D (xmin, xn, xmax) (ymin, yn, ymax) = DF.fillBuilder buildr
  where
    binX = DH.binD xmin xn xmax
    binY = DH.binD ymin yn ymax
    bins = Bin2D.Bin2D binX binY
    buildr = DF.mkSimple bins

histo :: (Foldable v, Unbox a, Num a) =>
         (Double, Int, Double)
      -> v Double
      -> DH.Histogram DF.BinD a
histo (xmin, n, xmax) = DF.fillBuilder buildr
  where
    bins = DH.binD xmin n xmax
    buildr = DF.mkSimple bins

-- A SIMPLE NEURAL NETWORK
-- Model Setup

data NN
  = NN
      { biass :: Vector Double,
        weights :: Vector Double,
        sigma :: Double
      }
  deriving (Eq, Show)

data Data
  = Data
      { xValue :: Double,
        yValue :: Double
      }
  deriving (Eq, Show)

forwardNN :: NN -> Double -> Double
forwardNN (NN bs ws _) x =
    ws `dot` cmap activation (scalar x - bs)
  where activation x = if x < 0 then 0 else 1

errorModel :: Double -> Double -> Double -> Log Double
errorModel mean std = normalPdf mean std

likelihood :: NN -> Data -> Log Double
likelihood nn (Data xObs yObs) =
    errorModel yMean ySigma yObs
  where
    ySigma = sigma nn
    yMean = forwardNN nn xObs

-- Prior, Posterior and Predictive Distribution

postNN :: MonadInfer m => m NN -> [Data] ->  m NN
postNN pr obs = do
  nn <- pr
  forM_ obs (score . likelihood nn)
  return nn

uniformVec :: MonadSample m => (Double, Double) -> Int -> m (Vector Double)
uniformVec (wmin, wmax) nelements =
  vector <$> replicateM nelements (uniform wmin wmax)

priorNN :: MonadSample m => Int -> m NN
priorNN nnodes = do
  bias <- uniformVec (0, 10) nnodes
  weight <- uniformVec (-10, 10) nnodes
  sigma <- uniform 0.5 1.5
  return $ NN bias weight sigma

predDist :: MonadInfer m => m NN -> m (NN, Data)
predDist pr = do
  nn <- pr
  x <- uniform 0 10
  y <- uniform (-5) 10
  score $ likelihood nn (Data x y)
  return (nn, Data x y)

max_stars_repo_name: leduin/Bayesian-Neural-Network

max_issues_repo_name: leduin/Bayesian-Neural-Network

max_forks_repo_name: leduin/Bayesian-Neural-Network


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

module Spiral.RootOfUnity (
    RootOfUnity(..)
  ) where

import Data.Complex
import Data.Complex.Cyclotomic

-- | Types that have primitive roots of unity.
class Fractional a => RootOfUnity a where
    -- | @rootOfUnity n k@ computes the primitive @n@th root of unity raised to
    -- the @k@th power.
    rootOfUnity :: Int -> Int -> a
    rootOfUnity n k = omega n ^^ k

    -- | Computes the primitive @n@th root of unity.
    omega :: Int -> a
    omega n = rootOfUnity n 1

instance RealFloat a => RootOfUnity (Complex a) where
    -- | $e^{-2 \pi i \frac{k}{n}$
    rootOfUnity n k = exp (-2*pi*i/fromIntegral n)^^k
      where
        i = 0:+1

instance RootOfUnity Cyclotomic where
    -- 'Data.Complex.Cyclotomic.e' computes $e^{2 \pi i/n}$, so we need to take
    -- the reciprocal.
    rootOfUnity n k = (1/e (fromIntegral n))^^k

max_stars_repo_name: mainland/hspiral

max_issues_repo_name: mainland/hspiral

max_forks_repo_name: mainland/hspiral


text:
{-# LANGUAGE DeriveDataTypeable        #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE FlexibleInstances         #-}
{-# LANGUAGE FunctionalDependencies    #-}
{-# LANGUAGE MultiParamTypeClasses     #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE RecordWildCards           #-}
{-# LANGUAGE TypeFamilies              #-}
{-# LANGUAGE UndecidableInstances      #-}

-----------------------------------------------------------------------------
-- |
-- Module      :  Plots.Types.Histogram
-- Copyright   :  (C) 2015 Christopher Chalmers
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Christopher Chalmers
-- Stability   :  experimental
-- Portability :  non-portable

-- A histogram is a graphical representation of the distribution of
-- numerical data. It is an estimate of the probability distribution of
-- a continuous variable.
--
----------------------------------------------------------------------------

module Plots.Types.Histogram
  (
    -- * Histogram plot
    HistogramPlot

    -- ** Already computed histograms
  , computedHistogram

    -- ** Histogram options
  , HistogramOptions
  , HasHistogramOptions (..)

    -- ** Normalisation
  , NormalisationMethod
  , count
  , probability
  , countDensity
  , pdf
  , cumilative
  , cdf

    -- ** Plotting histograms
  , histogramPlot
  , histogramPlot'
  , histogramPlotOf
  , histogramPlotOf'

    -- * Low level constructors
  , mkComputedHistogram
  , mkHistogramPlot
  ) where

import           Control.Monad.State.Lazy

import qualified Data.Foldable               as F
import           Data.Maybe
import           Data.Typeable

import qualified Data.Vector                 as V
import qualified Statistics.Sample.Histogram as Stat

import           Diagrams.Core.Transform     (fromSymmetric)
import           Diagrams.Prelude
import           Linear.V2                   (_yx)

import           Plots.Axis
import           Plots.Style
import           Plots.Types
import           Plots.Util

-- | Construct a rectangle of size $v$ with the bottom left at point $p$.
rectBL :: (InSpace V2 n t, TrailLike t) => Point V2 n -> V2 n -> t
rectBL p (V2 x y) =
  trailLike $ fromOffsets [V2 x 0, V2 0 y, V2 (-x) 0] # closeTrail `at` p

------------------------------------------------------------------------
-- GHistogram plot
------------------------------------------------------------------------

-- | Simple histogram type supporting uniform bins.
data HistogramPlot n = HistogramPlot
  { hWidth  :: n
  , hStart  :: n
  , hValues :: [n]
  , hOrient :: Orientation
  } deriving Typeable

type instance V (HistogramPlot n) = V2
type instance N (HistogramPlot n) = n

instance OrderedField n => Enveloped (HistogramPlot n) where
  getEnvelope HistogramPlot {..} =
    -- don't like this redundant code
    getEnvelope . orient hOrient _reflectXY id . (id :: Path v n -> Path v n) $
      ifoldMap drawBar hValues
    where
      drawBar i h = rectBL (mkP2 x 0) (V2 hWidth h)
        where x = hStart + fromIntegral i * hWidth

instance (TypeableFloat n, Renderable (Path V2 n) b)
    => Plotable (HistogramPlot n) b where
  renderPlotable s sty HistogramPlot {..} =
    ifoldMap drawBar hValues
      # orient hOrient _reflectXY id
      # applyAreaStyle sty
      # transform (s^.specTrans)
    where
      drawBar i h = rectBL (mkP2 x 0) (V2 hWidth h)
        where x = hStart + fromIntegral i * hWidth

  defLegendPic sty HistogramPlot {..}
    = centerXY
    . applyAreaStyle sty'
    . orient hOrient _reflectXY id
    $ alignB (rect 4 7) ||| alignB (rect 4 10) ||| alignB (rect 4 6)
    where
      -- The legend bars don't look right if the line width is too big so we limit it
      sty' = sty & areaStyle . _lw %~ atMost (local 0.8)

instance HasOrientation (HistogramPlot n) where
  orientation = lens hOrient $ \hp o -> hp {hOrient = o}

------------------------------------------------------------------------
-- Simple histogram plot
------------------------------------------------------------------------

-- | Plot an already computed histogram with equally sized bins.
computedHistogram
  :: (MonadState (Axis b V2 n) m,
      Plotable (HistogramPlot n) b,
      F.Foldable f)
  => n   -- ^ start of first bin
  -> n   -- ^ width of each bin
  -> f n -- ^ heights of the bins
  -> State (Plot (HistogramPlot n) b) ()
  -> m ()
computedHistogram x0 w xs = addPlotable (mkComputedHistogram x0 w xs)

-- | Construct a 'HistogramPlot' from raw histogram data.
mkComputedHistogram
  :: F.Foldable f
  => n -- ^ start of first bin
  -> n -- ^ width of each bin
  -> f n -- ^ heights of the bins
  -> HistogramPlot n
mkComputedHistogram x0 w xs = HistogramPlot x0 w (F.toList xs) Horizontal

----------------------------------------------------------------------------
-- Building histograms
----------------------------------------------------------------------------

-- example setup
-- > import Plots
-- > sampleData :: [Double]
-- > sampleData =
-- >   [5.1,4.9,4.7,4.6,5.0,5.4,4.6,5.0,4.4,4.9
-- >   ,5.4,4.8,4.8,4.3,5.8,5.7,5.4,5.1,5.7,5.1
-- >   ,5.4,5.1,4.6,5.1,4.8,5.0,5.0,5.2,5.2,4.7
-- >   ,4.8,5.4,5.2,5.5,4.9,5.0,5.5,4.9,4.4,5.1
-- >   ,5.0,4.5,4.4,5.0,5.1,4.8,5.1,4.6,5.3,5.0
-- >   ,7.0,6.4,6.9,5.5,6.5,5.7,6.3,4.9,6.6,5.2
-- >   ,5.0,5.9,6.0,6.1,5.6,6.7,5.6,5.8,6.2,5.6
-- >   ,5.9,6.1,6.3,6.1,6.4,6.6,6.8,6.7,6.0,5.7
-- >   ,5.5,5.5,5.8,6.0,5.4,6.0,6.7,6.3,5.6,5.5
-- >   ,5.5,6.1,5.8,5.0,5.6,5.7,5.7,6.2,5.1,5.7
-- >   ,6.3,5.8,7.1,6.3,6.5,7.6,4.9,7.3,6.7,7.2
-- >   ,6.5,6.4,6.8,5.7,5.8,6.4,6.5,7.7,7.7,6.0
-- >   ,6.9,5.6,7.7,6.3,6.7,7.2,6.2,6.1,6.4,7.2
-- >   ,7.4,7.9,6.4,6.3,6.1,7.7,6.3,6.4,6.0,6.9
-- >   ,6.7,6.9,5.8,6.8,6.7,6.7,6.3,6.5,6.2,5.9
-- >   ]
--
-- > mkNmExample nm = r2Axis &~ do
-- >   yMin ?= 0
-- >   histogramPlot sampleData $ do
-- >     normaliseSample .= nm
-- > countDia = renderAxis $ mkNmExample count
-- > probabilityDia = renderAxis $ mkNmExample probability
-- > countDensityDia = renderAxis $ mkNmExample countDensity
-- > pdfDia = renderAxis $ mkNmExample pdf
-- > cumilativeDia = renderAxis $ mkNmExample cumilative
-- > cdfDia = renderAxis $ mkNmExample cdf

-- Histogram options ---------------------------------------------------

-- | The way to normalise the data from a histogram. The default method
--   is 'count'.
newtype NormalisationMethod =
  NM { runNM :: forall n. Fractional n => n -> V.Vector n -> V.Vector n }
 -- width -> heights -> normalised heights

instance Default NormalisationMethod where
  def = count

-- | The height of each bar is the number of observations. This is the
--   'Default' method.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_countDia.svg#diagram=countDia&height=350>>
count :: NormalisationMethod
count = NM $ \_ v -> v

-- | The sum of the heights of the bars is equal to 1.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_probabilityDia.svg#diagram=probabilityDia&height=350>>
probability :: NormalisationMethod
probability = NM $ \_ v -> v ^/ V.sum v

-- | The height of each bar is @n / w@ where @n@ is the number of
--   observations and @w@ is the total width.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_countDensityDia.svg#diagram=countDensityDia&height=350>>
countDensity :: NormalisationMethod
countDensity = NM $ \w v -> v ^/ w

-- | The total area of the bars is @1@. This gives a probability density
--   function estimate.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_pdfDia.svg#diagram=pdfDia&height=350>>
pdf :: NormalisationMethod
pdf = NM $ \w v -> v ^/ (w * V.sum v)

-- | The height of each bar is the cumulative number of observations in
--   each bin and all previous bins. The height of the last bar is the
--   total number of observations.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_cumilativeDia.svg#diagram=cumilativeDia&height=350>>
cumilative :: NormalisationMethod
cumilative = NM $ \_ -> V.scanl1 (+)

-- | Cumulative density function estimate. The height of each bar is
--   equal to the cumulative relative number of observations in the bin
--   and all previous bins. The height of the last bar is 1.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_cdfDia.svg#diagram=cdfDia&height=350>>
cdf :: NormalisationMethod
cdf = NM $ \_ v -> V.scanl1 (+) v ^/ V.sum v

-- | Options for binning histogram data. For now only very basic
--   histograms building is supported.
data HistogramOptions n = HistogramOptions
  { hBins   :: Int
  , hRange  :: Maybe (n, n)
  , hNorm   :: NormalisationMethod
  , oOrient :: Orientation
  }

type instance V (HistogramOptions n) = V2
type instance N (HistogramOptions n) = n

instance Default (HistogramOptions n) where
  def = HistogramOptions
    { hBins   = 10
    , hRange  = Nothing
    , hNorm   = def
    , oOrient = Vertical
    }

instance HasOrientation (HistogramOptions n) where
  orientation = lens oOrient $ \ho o -> ho {oOrient = o}

class HasOrientation a => HasHistogramOptions a where
  -- | Options for building the histogram from data.
  histogramOptions :: Lens' a (HistogramOptions (N a))

  -- | The number of bins (bars) to use for the histogram. Must be
  --   positive.
  --
  --   'Default' is @10@.
  numBins :: Lens' a Int
  numBins = histogramOptions . lens hBins (\ho n -> ho {hBins = n})

  -- | The range of data to consider when building the histogram. Any
  --   data outside the range is ignored.
  --
  --   'Default' is 'Nothing'.
  binRange :: Lens' a (Maybe (N a, N a))
  binRange = histogramOptions . lens hRange (\ho r -> ho {hRange = r})

  -- | Should the resulting histogram be normalised so the total area is
  --   1.
  --
  --   'Default' is False.
  normaliseSample :: Lens' a NormalisationMethod
  normaliseSample = histogramOptions . lens hNorm (\ho b -> ho {hNorm = b})

instance HasHistogramOptions (HistogramOptions n) where
  histogramOptions = id

instance HasHistogramOptions a => HasHistogramOptions (Plot a b) where
  histogramOptions = rawPlot . histogramOptions

-- | Create a histogram by binning the data using the
--   'HistogramOptions'.
mkHistogramPlot
  :: (F.Foldable f, RealFrac n)
  => HistogramOptions n -> f n -> HistogramPlot n
mkHistogramPlot HistogramOptions {..} xs =
  HistogramPlot
    { hWidth  = w
    , hStart  = a
    , hValues = V.toList $ runNM hNorm w ns
    , hOrient = Vertical
    }
  where
    w     = (b - a) / fromIntegral hBins
    ns    = Stat.histogram_ hBins a b v
    v     = V.fromList (F.toList xs)
    (a,b) = fromMaybe (range hBins v) hRange

-- Taken from Statistics, which was limited to 'Double'.
range :: (Ord n, Fractional n)
      => Int                    -- ^ Number of bins (must be positive).
      -> V.Vector n             -- ^ Sample data (cannot be empty).
      -> (n, n)
range nBins xs
    | nBins < 1 = error "Plots.Types.Histogram: invalid bin count"
    | V.null xs = error "Plots.Types.Histogram: empty sample"
    | lo == hi  = case abs lo / 10 of
                    a | a < 1e-6   -> (-1,1)
                      | otherwise  -> (lo - a, lo + a)
    | otherwise = (lo-d, hi+d)
  where
    d | nBins == 1 = 0
      | otherwise  = (hi - lo) / ((fromIntegral nBins - 1) * 2)
    (lo,hi)        = minMaxOf folded xs
{-# INLINE range #-}

-- |
-- mkWeightedHistogram
--   :: (F.Foldable f, OrderdField n)
--   => HistogramOptions n -> [(n, n)] -> HistogramPlot n
-- mkWeightedHistogram

------------------------------------------------------------------------
-- Histogram
------------------------------------------------------------------------

-- $ histogram
-- Histograms display data as barplot of x data, bin y data.
-- Box plots have the following lenses:
--
-- @
-- * 'setBin' :: 'Lens'' ('BoxPlot' v n) 'Double' - 10
-- @

-- | Add a 'HistogramPlot' to the 'AxisState' from a data set.
--
-- === __Example__
--
-- <<diagrams/src_Plots_Types_Histogram_histogramExample.svg#diagram=histogramExample&height=350>>
--
-- > import Plots
-- > histogramAxis :: Axis B V2 Double
-- > histogramAxis = r2Axis &~ do
-- >   histogramPlot sampleData $ do
-- >     key "histogram"
-- >     plotColor .= blue
-- >     areaStyle . _opacity .= 0.5
--
-- > histogramExample = renderAxis histogramAxis
histogramPlot
  :: (MonadState (Axis b V2 n) m, Plotable (HistogramPlot n) b, F.Foldable f, RealFrac n)
  => f n -- ^ data
  -> State (Plot (HistogramOptions n) b) () -- ^ changes to plot options
  -> m () -- ^ add plot to axis
histogramPlot ns s = addPlot (hoPlot & rawPlot %~ \ho -> mkHistogramPlot ho ns)
  where hoPlot = mkPlot def &~ s

-- | Make a 'HistogramPlot' without changes to the plot options.
histogramPlot'
  :: (MonadState (Axis b V2 n) m, Plotable (HistogramPlot n) b, F.Foldable f, RealFrac n)
  => f n -- ^ data
  -> m () -- ^ add plot to axis
histogramPlot' d = histogramPlot d (return ())

-- | Add a 'HistogramPlot' using a fold over the data.
histogramPlotOf
  :: (MonadState (Axis b V2 n) m, Plotable (HistogramPlot n) b, RealFrac n)
  => Fold s n -- ^ fold over the data
  -> s        -- ^ data to fold
  -> State (Plot (HistogramOptions n) b) () -- ^ change to the plot
  -> m () -- ^ add plot to the 'Axis'
histogramPlotOf f s = histogramPlot (toListOf f s)

-- | Same as 'histogramPlotOf' without any changes to the plot.
histogramPlotOf'
  :: (MonadState (Axis b V2 n) m, Plotable (HistogramPlot n) b, RealFrac n)
  => Fold s n -> s -> m ()
histogramPlotOf' f s = histogramPlotOf f s (return ())

-- temporary functions that will be in next lib release

_reflectionXY :: (Additive v, R2 v, Num n) => Transformation v n
_reflectionXY = fromSymmetric $ (_xy %~ view _yx) <-> (_xy %~ view _yx)

_reflectXY :: (InSpace v n t, R2 v, Transformable t) => t -> t
_reflectXY = transform _reflectionXY

max_stars_repo_name: adamConnerSax/plots

max_issues_repo_name: adamConnerSax/plots

max_forks_repo_name: adamConnerSax/plots


text:
module Raycast (intersect, findHits, raycast) where

import           Control.Monad (guard)
import           Data.Maybe (isJust, fromJust, listToMaybe)
import Types
import Ray (createHit)
import Numeric.Vector
import Numeric.Scalar

-- https://courses.cs.washington.edu/courses/csep557/10au/lectures/triangle_intersection.pdf
intersect :: Ray -> Object -> Tri -> Maybe RayHit
intersect ray@(Ray pos dir) obj tri = do
  let a = v_getPos $ t_getP0 tri
      b = v_getPos $ t_getP1 tri
      c = v_getPos $ t_getP2 tri
      n = normalized $ (b - a) × (c - a)
      d = n · a
      t = (d - n · pos) / (n · dir)
      q = pos + dir * (fromScalar t)

      areaQAB = ((b - a) × (q - a)) · n
      areaQBC = ((c - b) × (q - b)) · n
      areaQCA = ((a - c) × (q - c)) · n
      areaABC = ((b - a) × (c - a)) · n

  guard $ areaQAB >= 0
  guard $ areaQBC >= 0
  guard $ areaQCA >= 0
  guard $ t > 0
  -- Only counts as a hit if the ray is traveling towards the front of the triangle
  -- (when direction and normal are closer to antiparallel)
  guard $ (dir · n) < 0

  let alpha = unScalar $ areaQBC / areaABC
      beta = unScalar $ areaQCA / areaABC
      gamma = unScalar $ areaQAB / areaABC

  let hitPos = pos + dir * (fromScalar t)

  return $ createHit obj tri ray (unScalar t) (alpha, beta, gamma) hitPos

findHits :: Scene -> Ray -> [RayHit]
findHits s ray = let objs = s_getObjects s
                     intersections = concatMap (\o -> map (intersect ray o) $ o_getTris o) objs
                 -- It's safe to use `fromJust` here because the filter guarantees all values are `Just`
                 in  map fromJust $ filter isJust intersections

minimumBy :: (Ord b) => (a -> b) -> [a] -> a
minimumBy f xs = foldr1 (\m x -> if f x < f m then x else m) xs

raycast :: Scene -> Ray -> Maybe RayHit
raycast s ray@(Ray origin _) = let hs = findHits s ray
                               in  if   length hs == 0 then Nothing
                                   else Just $ minimumBy rh_getDistance hs


text:
module LatticeSymmetries.IO where

import Data.Aeson
import Data.Aeson.Types (typeMismatch)
import Data.Complex
import Data.Scientific

data SymmetrySpec = SymmetrySpec !(NonEmpty Int) !Int
  deriving stock (Read, Show, Eq)

instance FromJSON SymmetrySpec where
  parseJSON = withObject "symmetry" $ \v ->
    SymmetrySpec
      <$> v .: "permutation"
      <*> v .: "sector"

instance ToJSON SymmetrySpec where
  toJSON (SymmetrySpec permutation sector) =
    object ["permutation" .= permutation, "sector" .= sector]

data BasisSpec = BasisSpec !Int !(Maybe Int) !(Maybe Int) ![SymmetrySpec]
  deriving stock (Read, Show, Eq)

instance FromJSON BasisSpec where
  parseJSON = withObject "basis" $ \v ->
    BasisSpec
      <$> v .: "number_spins"
      <*> v .:? "hamming_weight"
      <*> v .:? "spin_inversion"
      <*> v .:! "symmetries" .!= []

instance ToJSON BasisSpec where
  toJSON (BasisSpec numberSpins hammingWeight spinInversion symmetries) =
    object
      [ "number_spins" .= numberSpins,
        "hamming_weight" .= maybe Null toJSON hammingWeight,
        "spin_inversion" .= maybe Null toJSON spinInversion,
        "symmetries" .= symmetries
      ]

newtype WrappedBasisSpec = WrappedBasisSpec BasisSpec

instance FromJSON WrappedBasisSpec where
  parseJSON = withObject "config" $ \v -> WrappedBasisSpec <$> v .: "basis"

instance FromJSON (Complex Double) where
  parseJSON (Number x) = pure . fromReal . toRealFloat $ x
    where
      fromReal :: Num a => a -> Complex a
      fromReal x' = x' :+ 0
  parseJSON v@(Array xs) = case (toList xs) of
    [re, im] -> (:+) <$> parseJSON re <*> parseJSON im
    _ -> typeMismatch "Complex" v
  parseJSON v = typeMismatch "Complex" v

newtype DenseMatrixSpec = DenseMatrixSpec [[Complex Double]]
  deriving stock (Read, Show, Eq)
  deriving newtype (FromJSON)

data InteractionSpec = InteractionSpec ![[Complex Double]] ![[Int]]
  deriving stock (Read, Show, Eq)

instance FromJSON InteractionSpec where
  parseJSON = withObject "interaction" $ \v ->
    InteractionSpec
      <$> v .: "matrix"
      <*> v .: "sites"

data OperatorSpec = OperatorSpec !Text !(NonEmpty InteractionSpec)
  deriving stock (Read, Show)

instance FromJSON OperatorSpec where
  parseJSON = withObject "operator" $ \v ->
    OperatorSpec
      <$> v .: "name"
      <*> v .: "terms"

data ConfigSpec = ConfigSpec !BasisSpec !OperatorSpec
  deriving stock (Read, Show)

instance FromJSON ConfigSpec where
  parseJSON = withObject "config" $ \v ->
    ConfigSpec
      <$> v .: "basis"
      <*> v .: "hamiltonian"

max_stars_repo_name: twesterhout/lattice-symmetries-haskell

max_issues_repo_name: twesterhout/lattice-symmetries-haskell

max_forks_repo_name: twesterhout/lattice-symmetries-haskell


text:
{-# LANGUAGE BangPatterns #-}
module STC where

import           Control.Monad.Parallel      as MP
import           Data.Array.IArray           as IA
import           Data.Array.Repa             as R
import           Data.Complex
import           Data.List                   as L
import           Data.Vector.Storable        as VS
import           Data.Vector.Unboxed         as VU
import           FokkerPlanck.FourierSeries
import           FokkerPlanck.GreensFunction
import           FokkerPlanck.Histogram
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Parallel
import           Utils.Time

main = do
  args@(deviceIDsStr:numPointStr:deltaStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:deltaLogStr:taoStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:stdStr:numBatchR2FreqsStr:numBatchOriStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoint = read numPointStr :: Int
      delta = read deltaStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tao = read taoStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      numThread = read numThreadStr :: Int
      folderPath = "output/test/FourierSeriesOfPinwheels"
      maxScale = read maxScaleStr :: Double
      halfLogPeriod = log maxScale
      deltaLog = read deltaLogStr :: Double
      std = read stdStr :: Double
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      numBatchOri = read numBatchOriStr :: Int
      radius = fromIntegral $ div numPoint 2
  removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  printCurrentTime "Start computing coefficients..."
  hist <-
    computeFourierCoefficients
      deviceIDs
      numPoints
      numOrientation
      delta
      initScale
      thetaSigma
      tao
      numR2Freq
      periodR2
      std
      numBatchR2Freqs
      numBatchOri
      phiFreq
      rhoFreq
      thetaFreq
      scaleFreq
  hist <-
    sampleCartesian
      folderPath
      histFilePath
      deviceIDs
      radius
      144
      0.25 -- deltaLog
      initScale
      thetaSigma
      tao
      phiFreqs
      rhoFreqs
      thetaFreqs
      scaleFreqs
  printCurrentTime "Done"
  printCurrentTime "Start Convloution.."
  let !coefficients =
        normalizeFreqArr' std phiFreqs rhoFreqs . getNormalizedHistogramArr $
        hist
      initialDistribution =
        computeInitialDistributionFull' r2Freq radius phiFreq rhoFreq initPoints
      !harmonicsArray' =
        computeFourierSeriesOfLogPolarHarmonicsArray
          radius
          deltaLog
          r2Freq
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          halfLogPeriod :: IA.Array (Int, Int) (VS.Vector (Complex Double))
      -- !harmonicsArray' = computeHarmonicsArray numPoint 1 numPoint 1 phiFreqs rhoFreqs thetaFreqs scaleFreqs 1 radius
      source = convolveFull' coefficients harmonicsArray' initialDistribution
      inverseR2Harmonics =
        computeRectangularInverseHarmonics numPoint numPoint delta radius r2Freq :: [VS.Vector (Complex Double)]
      -- pinwheelVec =
      --   parMap
      --     rdeepseq
      --     (\(idx, vec) ->
      --        ( idx
      --        , VU.fromList . L.map (VS.sum . VS.zipWith (*) vec) $
      --          inverseR2Harmonics)) .
      --   assocs $
      --   harmonicsArray'
  -- print initialDistribution
  createDirectoryIfMissing True (folderPath </> "Pinwheels")
  plotImageRepaComplex (folderPath </> "Pinwheels" </> "Pinwheel_(1,1).png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numPoint :. numPoint) .
    VU.fromList . L.map (VS.sum . VS.zipWith (*) (harmonicsArray' IA.! (1, 1))) $
    inverseR2Harmonics
  -- MP.mapM_
  --   (\((a, b), vec) ->
  --      plotImageRepaComplex
  --        (folderPath </> "Pinwheels" </> printf "Pinwheel_(%d,%d).png" a b) .
  --      ImageRepa 8 . fromUnboxed (Z :. (1 :: Int) :. numPoint :. numPoint) $
  --      vec)
  --   pinwheelVec
  -- printCurrentTime "Done"
  printCurrentTime "Start ploting.."
  plotDFTArrayPower
    (folderPath </> "SourcePower.png")
    numPoint
    numPoint
    -- (2 * r2Freq + 1)
    -- (2 * r2Freq + 1)
    source
  plotFourierCoefficientsPower
    (folderPath </> "Source.png")
    numPoint
    numPoint
    inverseR2Harmonics
    source
  printCurrentTime "Done"

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 FunctionalDependencies #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE OverloadedLists #-}

module Data.Network
  (Network, Layer, Batch, Data, Layers(Output, (:=>)),
   makeRandomLayer, makeRandomNetwork, run, compute, runLayers, runNetwork, fromList, toList
  )
  --   makeRandomNetwork, makeRandomModel, makeModel, fromList,
  --   runLayer, runLayer', runNetwork, runLayers, runLayers',
  --   runLayerBatch, runLayerBatch', runNetworkBatch, runLayersBatch, runLayersBatch',
  --   Input, Output, Data, Batch, BatchIn, BatchOut, Sample, BatchSample, Weights)
where

import GHC.Exts hiding ((<#))
import Prelude hiding (scanl)
import Data.Functor
-- import Data.Word8
import Data.Foldable hiding (toList)
import Data.Monoid
import System.Random
import Control.Monad.Random as R hiding (fromList)
import Control.Monad
import Control.DeepSeq
import Numeric.LinearAlgebra as LA hiding ((<>), fromList, toList)
import Numeric.LinearAlgebra.Data as LAD hiding (fromList, toList, size, R)
-- import Numeric.LinearAlgebra.Static as LAS hiding ((<>), fromList)
import Numeric.LinearAlgebra.Devel
import Data.Network.Utils (scanl)
import Data.Network.Activation

-- type ActivationFunction = forall a. Floating a => a -> a

type Input = Vector Double
type Output = Vector Double
type Data = Vector Double
type Weights = Matrix Double
type Batch = Matrix Double
type BatchIn = Matrix Double
type BatchOut = Matrix Double
type Sample = (Data, Data) -- Labeled data for supervised learning
type BatchSample = (Batch, Batch)
-- type Model a = (Network, ActivationFunction a) -- A network's structure with its activation function

-- Network layer = (weight matrix, bias vector)
-- Parameter a: type of data(e.g. Double)

-- What's a layer?
class (NFData a) => Layer i a | a -> i where
  compute :: a -> i -> i
  run     :: a -> i -> (i, i)

  makeRandomLayer :: (RandomGen g) => g -> ([Double] -> [Double]) -> Int -> Int -> ActivationFunction i -> a

-- What's a network?
class (NFData a) => Network i a | a -> i where
  runNetwork :: a -> i -> (i,i)
  runLayers  :: a -> i -> [(i,i)]
  makeRandomNetwork :: (RandomGen g) =>
    g -> ([Double] -> [Double]) -> (a -> a) -> Int -> [(Int, ActivationFunction i)] -> (Int, ActivationFunction i) -> a

-- Linked list Structure for feedforward neural network
data Layers :: * -> * where
  Output :: !a -> Layers a
  (:=>) :: !a -> !(Layers a) -> Layers a deriving (Eq)

infixr 5 :=>

instance (NFData a) => NFData (Layers a) where
  rnf (Output x) = rnf x
  rnf (l :=> ll) = rnf l `seq` rnf ll

instance Show a => Show (Layers a) where
  show (Output l) = "(Output) " ++ show l
  show (a :=> b) = show a ++ " -> " ++ show b

instance Functor Layers where
  fmap f (Output x) = (Output $ f x)
  fmap f (a :=> b) = (f a) :=> (fmap f b)

instance Foldable Layers where
  foldr f b (Output x) = f x b
  foldr f c (a :=> b) = f a (foldr f c b)

  foldl f a (Output x) = f a x
  foldl f a (b :=> c) = foldl f (f a b) c

instance Num a => Num (Layers a) where
  (Output l) + (Output r) = Output (l + r)
  (l :=> n) + (r :=> n') = (l + r) :=> (n + n')
  (Output l) - (Output r) = Output (l - r)
  (l :=> n) - (r :=> n') = (l - r) :=> (n - n')
  (Output l) * (Output r) = Output (l * r)
  (l :=> n) * (r :=> n') = (l * r) :=> (n * n')
  abs (Output l) = Output (abs l)
  abs (l :=> n) = (abs l :=> abs n)
  signum (Output l) = Output (signum l)
  signum (l :=> n) = (signum l) :=> (signum n)
  fromInteger a = Output (fromInteger a)

instance Fractional a => Fractional (Layers a) where
  (Output l) / (Output r) = Output (l/r)
  (l :=> n) / (r :=> n') = (l/r) :=> (n/n')
  recip (Output l) = Output (recip l)
  recip (l :=> n) = (recip l) :=> (recip n)
  fromRational a = Output (fromRational a)

    -- g -> ([Double] -> [Double]) -> Int -> [(Int, ActivationFunction i)] -> (Int, ActivationFunction i) -> a

instance IsList (Layers i) where
  type Item (Layers i) = i
  fromList []     = undefined
  fromList [x]    = Output x
  fromList (x:xs) = x :=> (fromList xs)

  toList (Output x) = [x]
  toList (x :=> y)  = x:(toList y)

-- fromList :: [a] -> Layers a
-- fromList [] = undefined
-- fromList [x] = Output x
-- fromList (x:xs) = x :=> (fromList xs)

-- make a fully connected layer; m = size of input(previous layer), n = size of output(next layer)
-- makeLayerFC :: Floating a => Int -> Int -> ActivationFunction a -> Layer a Double
-- makeLayerFC m n af = Layer { weights = (m><n) (repeat 0), biases = n |> (repeat 0), activation = af }

-- make a fully connected layer with randomly generated (normally distributed) weights and biases
-- makeLayerFCR :: MonadRandom m => RandDist -> Int -> Int -> ActivationFunction a -> m (Layer a Double)
-- makeLayerFCR dist m n = do
--   (seeds :: [Int]) <- getRandoms
--   (seed2 :: Int) <- getRandom
--   let w = fromRows $ fmap (\s -> randomVector s dist n) (take m seeds)
--       b = randomVector seed2 dist n
--   return $ Layer { weights = w, biases = b }

-- makeRandomNetwork :: MonadRandom m => RandDist -> Int -> [Int] -> Int -> m Network
-- makeRandomNetwork dist i [] o = Output <$> (makeLayerFCR dist i o)
-- makeRandomNetwork dist i (h:hs) o = (:=>) <$> makeLayerFCR dist i h <*> makeRandomNetwork dist h hs o

-- makeModel :: Network -> [ActivationFunction a] -> Model a
-- makeModel (Output l) (a:afs) = Output (l,a)
-- makeModel (l :=> n) (a:afs) = ((l,a) :=> makeModel n afs)

-- makeRandomModel :: MonadRandom m => RandDist -> Int -> [(Int, ActivationFunction a)] -> (Int, ActivationFunction a) -> m (Network i)
-- makeRandomModel dist i [] (o, oaf) = (Output . (,oaf)) <$> (makeLayerFCR dist i o)
-- makeRandomModel dist i ((h,af):hs) (o, oaf) = ((:=>) . (,af)) <$> makeLayerFCR dist i h <*> makeRandomModel dist h hs (o, oaf)

-- compute layer's net result (before activation/squashing)
-- computeLayer :: Layer a Double -> a -> Vector Double
-- computeLayer (Layer weights biases) input = biases + input <# weights

-- net linear combiner for a batch of input as a matrix, with rows as input vectors
-- computeLayerBatch :: Layer Batch Double -> Batch -> Batch
-- computeLayerBatch (Layer weights biases) inputs = (fromRows $ replicate (rows inputs) biases) + (inputs <> weights)

-- run layer with input vector
-- runLayer :: ActivationFunction Output -> Layer Double -> Input -> Output
-- runLayer af l i = activate af $ computeLayer l i

-- run layer with output vector, outputing both result and derivative
-- runLayer' :: ActivationFunction Output -> Layer Double -> Input -> (Output, Output)
-- runLayer' af l i = unAct af $ computeLayer l i

-- Batch variants

-- runLayerBatch :: ActivationFunction Batch -> Layer Double -> BatchIn -> BatchOut
-- runLayerBatch af l i = activate af $ computeLayerBatch l i

-- runLayerBatch' :: ActivationFunction Batch -> Layer Double -> BatchIn -> (BatchOut, BatchOut)
-- runLayerBatch' af l i = unAct af $ computeLayerBatch l i

-- run network on input vector with an activation function
-- runNetwork :: Model Data -> Input -> Output
-- runNetwork n i = foldl' ((flip . uncurry . flip) runLayer) i n

-- runNetworkBatch :: Model Batch -> BatchIn -> BatchOut
-- runNetworkBatch n i = foldl' ((flip . uncurry . flip) runLayerBatch) i n

-- run network on input, accumulating each layer's intermediary result,
-- with the first element of the resulting list being the last layer's output, ...
-- runLayers :: Model Data -> Input -> [Output]
-- runLayers n i = foldl' (\os@(o:os') (l,af) -> (runLayer af l o):os) [i] n

-- Same as above but with derivatives
-- runLayers' :: Model Data -> Input -> [(Output, Output)]
-- runLayers' n i = foldl' (\os@((o, o'):os') (l, af) -> (runLayer' af l o):os) [(i, mempty)] n

-- Batch variants
-- runLayersBatch :: Model Batch -> BatchIn -> [BatchOut]
-- runLayersBatch n i = foldl' (\os@(o:os') (l, af) -> (runLayerBatch af l o):os) [i] n

-- runLayersBatch' :: Model Batch -> BatchIn -> [(BatchOut, BatchOut)]
-- runLayersBatch' n i = foldl' (\os@((o, o'):os') (l, af) -> (runLayerBatch' af l o):os) [(i, mempty)] n

max_stars_repo_name: DrPyser/NeuralNetwork

max_issues_repo_name: DrPyser/NeuralNetwork

max_forks_repo_name: DrPyser/NeuralNetwork


text:
module Pfe3Metrics where
import Data.Maybe(mapMaybe)

import PFE0(allModules,pput)
import PFE3(parseModule)
import HsDeclStruct(DI(..))
import HsModule(hsModDecls)
import DefinedNames(definedType)
import HasBaseStruct(basestruct)
import Statistics

import PfeParse
import MUtils

pfe3MetricsCmds =
  [("classmetrics", (noArgs classMetrics,"number of instances per class metrics"))]
classMetrics =
  do ms <- allModules
     (classes,insts) <- apBoth concat . unzip # mapM getClassInstDecls ms
     let cinstcnt = [(c,length [()|i<-insts,i==c])|c<-classes]
     pput (ppStatistics "number of instances" "class" cinstcnt)

getClassInstDecls m = cls_insts . hsModDecls . snd # parseModule m

cls_insts ds = (mapMaybe className ds,concatMap instClass ds)

className d =
  case basestruct d of
    Just (HsClassDecl _ _ tp _ _) -> Just (definedType tp)
    _ -> Nothing

instClass d =
  case basestruct d of
    Just (HsInstDecl _ _ _ tp _) -> [definedType tp]
    Just (HsDataDecl _ _ _ _ cls) -> cls
    Just (HsNewTypeDecl _ _ _ _ cls) -> cls
    _ -> []

max_stars_repo_name: bubba/HaRe

max_issues_repo_name: bubba/HaRe

max_forks_repo_name: bubba/HaRe


text:
module Main where

import qualified Numeric.LinearAlgebra.Data

import Optimization.Types
import Optimization.NelderMead

trivialOpt :: Vector Double -> Double
trivialOpt v
  | length v != 3 = error "Vector is not length-3!"
  | otherwise     = ((v ! 0)**2) * (((v ! 1) - 1)**2) * (((v ! 2) - 2) ** 2)

main :: IO ()
main = do
  let
    tc      =  terminateAtOrAfter 0.001 1000
    optimum =  nelderMead tc trivialOpt (fromList [10.0, 12.1, -4.4])

  putStrLn $ show optimum

max_stars_repo_name: agbrooks/buzzwords

max_issues_repo_name: agbrooks/buzzwords

max_forks_repo_name: agbrooks/buzzwords


text:
module Rebase.Data.Complex
(
  module Data.Complex
)
where

import Data.Complex

max_stars_repo_name: felixonmars/rebase

max_issues_repo_name: felixonmars/rebase

max_forks_repo_name: felixonmars/rebase


text:
{-# LANGUAGE QuasiQuotes, MultiParamTypeClasses #-}
{-# LANGUAGE TemplateHaskell #-}
module Lib where

import qualified Language.C.Inline as C
import Foreign.C.Types
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA
import qualified Numeric.LinearAlgebra.Data as D
import qualified Numeric.LinearAlgebra.Devel as D
import Control.Arrow ((***))
import Data.Tuple (swap)
C.context (C.baseCtx <> C.vecCtx)

C.include "<math.h>"
C.include "<osqp/osqp.h>"

testsp = D.mkSparse [((0,0), 3.0), ((1,1), 2.0)]
testcsr = D.mkCSR [((0,0), 3.0), ((1,1), 2.0)]

-- getCSR (D.SparseR csr ncols nrows) =

assocVec :: (Num a, Ord a) => [a] -> [(Int,a)]
assocVec v = filter (\(_,x) -> x /= 0) (zipWith (,) [0..] v)

-- conveniently, will assume that if your list stops, it is implied the rest of it is zero. Or maybe inconvenitly from a type safety persepctive
assocMat :: (Num a, Ord a) => [[a]] -> [((Int,Int),a)]
assocMat m = let m' = map assocVec m in addRows m' where addRows m'' = concat $ zipWith (\i xs -> map (\(j,x) -> ((i,j),x)) xs) [0 ..] m''
-- AssocMatrix = [((Int, Int), Double)]

transposeAssoc = map (swap *** id)

mkGCSC = D.tr . D.mkSparse . transposeAssoc
mkGCSC' = mkGCSC . assocMat
mkCSC :: D.AssocMatrix -> CSC'
mkCSC = transpose'' . D.mkCSR . transposeAssoc
mkCSC' :: [[Double]] -> CSC'
mkCSC' = mkCSC . assocMat
-- The CSC type is internal to hmatrix.
data CSC' = CSC'
        { cscVals  :: V.Vector Double
        , cscRows  :: V.Vector CInt
        , cscCols  :: V.Vector CInt
        , cscNRows :: Int
        , cscNCols :: Int
        } deriving Show

-- data OSQPCSC = CSC  -- {nzmax :: CInt, m :: CInt, n :: CInt, col_p :: V.Vector Ptr, col  }
{-
instance LA.Transposable D.CSR CSC'
  where
  tr (D.CSR vs cs rs n m) = CSC' vs cs rs m n
  tr' = tr

instance LA.Transposable CSC' D.CSR
  where
  tr (CSC' vs rs cs n m) = D.CSR vs rs cs m n
  tr' = tr
-}
transpose'' (D.CSR vs cs rs n m) = CSC' vs cs rs m n
transpose''' (CSC' vs rs cs n m) = D.CSR vs rs cs m n

funboy :: IO ()
funboy = do
  x <- [C.exp| double{ cos(1) } |]
  print x
someFunc :: IO ()
someFunc = putStrLn "someFunc"

myfun :: CDouble -> CDouble
myfun x = [C.pure| double{ cos( $(double x)) } |]

-- exampleOSQPData =

-- Actually CDoubles
--data OSQPData = OSQPData {dn :: CInt, dm :: CInt, p :: CSC ,
--     a :: CSC, q :: V.Vector CDouble, l :: V.Vector CDouble, u :: V.Vector CDouble}

-- OSQPWorkspace *osqp_setup(const OSQPData *data, OSQPSettings *settings)
-- c_int osqp_solve(OSQPWorkspace *work)
-- c_int osqp_cleanup(OSQPWorkspace *work)

-- Maybe the easiest thing is to just inject default settings and marshal only what we need
-- Get a demo code in full c and start pulling pieces out.
-- l and u are the easiest

demo :: IO CInt
demo = [C.block|
int {
    // Load problem data
    c_float P_x[4] = {4.00, 1.00, 1.00, 2.00, };
    c_int P_nnz = 4;
    c_int P_i[4] = {0, 1, 0, 1, };
    c_int P_p[3] = {0, 2, 4, };
    c_float q[2] = {1.00, 1.00, };
    c_float A_x[4] = {1.00, 1.00, 1.00, 1.00, };
    c_int A_nnz = 4;
    c_int A_i[4] = {0, 1, 0, 2, };
    c_int A_p[3] = {0, 2, 4, };
    c_float l[3] = {1.00, 0.00, 0.00, };
    c_float u[3] = {1.00, 0.70, 0.70, };
    c_int n = 2;
    c_int m = 3;

    // Problem settings
    OSQPSettings * settings = (OSQPSettings *)c_malloc(sizeof(OSQPSettings));

    // Structures
    OSQPWorkspace * work;  // Workspace
    OSQPData * data;  // OSQPData

    // Populate data
    data = (OSQPData *)c_malloc(sizeof(OSQPData));
    data->n = n;
    data->m = m;
    data->P = csc_matrix(data->n, data->n, P_nnz, P_x, P_i, P_p);
    data->q = q;
    data->A = csc_matrix(data->m, data->n, A_nnz, A_x, A_i, A_p);
    data->l = l;
    data->u = u;

    // Define Solver settings as default
    osqp_set_default_settings(settings);
    settings->alpha = 1.0; // Change alpha parameter

    // Setup workspace
    work = osqp_setup(data, settings);

    // Solve Problem
    osqp_solve(work);

    // Cleanup
    osqp_cleanup(work);
    c_free(data->A);
    c_free(data->P);
    c_free(data);
    c_free(settings);

    return 0;
} |]


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables   #-}
module Grenade.Layers.FullyConnected (
    FullyConnected (..)
  , FullyConnected' (..)
  , randomFullyConnected
  ) where

import           Control.Monad.Random hiding (fromList)

import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits

import qualified Numeric.LinearAlgebra as LA
import           Numeric.LinearAlgebra.Static

import           Grenade.Core

import           Grenade.Layers.Internal.Update

-- | A basic fully connected (or inner product) neural network layer.
data FullyConnected i o = FullyConnected
                        !(FullyConnected' i o)   -- Neuron weights
                        !(FullyConnected' i o)   -- Neuron momentum

data FullyConnected' i o = FullyConnected'
                         !(R o)   -- Bias
                         !(L o i) -- Activations

instance Show (FullyConnected i o) where
  show FullyConnected {} = "FullyConnected"

instance (KnownNat i, KnownNat o) => UpdateLayer (FullyConnected i o) where
  type Gradient (FullyConnected i o) = (FullyConnected' i o)

  runUpdate lp (FullyConnected (FullyConnected' oldBias oldActivations) (FullyConnected' oldBiasMomentum oldMomentum)) (FullyConnected' biasGradient activationGradient) =
    let (newBias, newBiasMomentum)    = descendVector (learningRate lp) (learningMomentum lp) (learningRegulariser lp) oldBias biasGradient oldBiasMomentum
        (newActivations, newMomentum) = descendMatrix (learningRate lp) (learningMomentum lp) (learningRegulariser lp) oldActivations activationGradient oldMomentum
    in FullyConnected (FullyConnected' newBias newActivations) (FullyConnected' newBiasMomentum newMomentum)

  createRandom = randomFullyConnected

instance (KnownNat i, KnownNat o) => Layer (FullyConnected i o) ('D1 i) ('D1 o) where
  type Tape (FullyConnected i o) ('D1 i) ('D1 o) = R i
  -- Do a matrix vector multiplication and return the result.
  runForwards (FullyConnected (FullyConnected' wB wN) _) (S1D v) = (v, S1D (wB + wN #> v))

  -- Run a backpropogation step for a full connected layer.
  runBackwards (FullyConnected (FullyConnected' _ wN) _) x (S1D dEdy) =
          let wB'  = dEdy
              mm'  = dEdy `outer` x
              -- calcluate derivatives for next step
              dWs  = tr wN #> dEdy
          in  (FullyConnected' wB' mm', S1D dWs)

instance (KnownNat i, KnownNat o) => Serialize (FullyConnected i o) where
  put (FullyConnected (FullyConnected' b w) _) = do
    putListOf put . LA.toList . extract $ b
    putListOf put . LA.toList . LA.flatten . extract $ w

  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy i)
      b     <- maybe (fail "Vector of incorrect size") return . create . LA.fromList =<< getListOf get
      k     <- maybe (fail "Vector of incorrect size") return . create . LA.reshape f . LA.fromList =<< getListOf get
      let bm = konst 0
      let mm = konst 0
      return $ FullyConnected (FullyConnected' b k) (FullyConnected' bm mm)

randomFullyConnected :: (MonadRandom m, KnownNat i, KnownNat o)
                     => m (FullyConnected i o)
randomFullyConnected = do
    s1    <- getRandom
    s2    <- getRandom
    let wB = randomVector  s1 Uniform * 2 - 1
        wN = uniformSample s2 (-1) 1
        bm = konst 0
        mm = konst 0
    return $ FullyConnected (FullyConnected' wB wN) (FullyConnected' bm mm)

max_stars_repo_name: LuisChDev/grenade

max_issues_repo_name: LuisChDev/grenade

max_forks_repo_name: LuisChDev/grenade


text:
module Reanimate.Morph.Triagulate where

-- import Graphics.SvgTree
import Linear.V2
import Linear.Vector
-- import Reanimate
import Numeric.LinearAlgebra
-- import Numeric.LinearAlgebra.Sparse
-- import Data.Sparse.SpMatrix
-- import Data.Sparse.SpVector

type Poly = [V2 Double]

-- renderPolygon :: Poly -> SVG
-- renderPolygon = undefined

triangulate :: Poly -> Poly -> ([Poly],[Poly])
triangulate = undefined

type Point = V2 Double

p1,p2,p3,p4 :: Point
p1 = V2 0 0
p2 = V2 2 0
p3 = V2 1 2
p4 = V2 1 1

-- let (t1,t2,t3) = barycentricCoordinates p1 p2 p3 p4
-- p4 = t1*p1 + t2*p2 + t3*p3
--
--
--                       [ t1
--                       , t2
--                       , t3 ]
-- [ [p1_x, p2_x, p3_x ] [ p4_x
-- , [p1_y, p2_y, p3_y ] , p4_y
-- , [1, 1, 1] ]         , 1 ]
-- barycentricCoordinates :: Point -> Point -> Point -> Point -> (Double, Double, Double)
barycentricCoordinates (V2 p1x p1y) (V2 p2x p2y) (V2 p3x p3y) (V2 p4x p4y) =
    linearSolve m v
    -- (m, v)
  where
    m = (3><3) [p1x, p2x, p3x, p1y, p2y, p3y, 1, 1, 1 :: Double]
    v = (3><1) [ p4x, p4y, 1 :: Double ]

bCoords (V2 x1 y1) (V2 x2 y2) (V2 x3 y3) (V2 x y) =
    (lam1, lam2, lam3)
  where
    lam1 = ((y2-y3)*(x-x3) + (x3 - x2)*(y-y3)) /
           ((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
    lam2 = ((y3-y1)*(x-x3) + (x1-x3)*(y-y3)) /
           ((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
    lam3 = 1 - lam1 - lam2

max_stars_repo_name: sureyeaah/reanimate

max_issues_repo_name: sureyeaah/reanimate

max_forks_repo_name: sureyeaah/reanimate


text:
{-# LANGUAGE TemplateHaskell #-}
{-|
    Module      :  Numeric.MixedType.Complex
    Description :  Instances for Data.Complex
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

    Maintainer  :  mikkonecny@gmail.com
    Stability   :  experimental
    Portability :  portable

    Instances for "Data.Complex".
-}

{-# OPTIONS_GHC -Wno-orphans #-}
module Numeric.MixedTypes.Complex
(
  tComplex
)
where

import Numeric.MixedTypes.PreludeHiding
-- import qualified Prelude as P
-- import Text.Printf

import Utils.TH.DeclForTypes

import Data.Complex

import Numeric.MixedTypes.Literals
import Numeric.MixedTypes.Bool
import Numeric.MixedTypes.Eq
import Numeric.MixedTypes.MinMaxAbs
import Numeric.MixedTypes.AddSub
import Numeric.MixedTypes.Mul
import Numeric.MixedTypes.Div
-- import Numeric.MixedTypes.Power
-- import Numeric.MixedTypes.Field
import Numeric.MixedTypes.Elementary

tComplex :: T t -> T (Complex t)
tComplex (T tName) = T ("(Complex " ++ tName ++ ")")

instance (ConvertibleExactly Integer t) => (ConvertibleExactly Integer (Complex t))
  where
  safeConvertExactly n =
    do
    nT <- safeConvertExactly n
    zT <- safeConvertExactly 0
    return $ nT :+ zT

instance (ConvertibleExactly Int t) => (ConvertibleExactly Int (Complex t))
  where
  safeConvertExactly n =
    do
    nT <- safeConvertExactly n
    zT <- safeConvertExactly (int 0)
    return $ nT :+ zT

instance (ConvertibleExactly Rational t) => (ConvertibleExactly Rational (Complex t))
  where
  safeConvertExactly r =
    do
    rT <- safeConvertExactly r
    zT <- safeConvertExactly 0.0
    return $ rT :+ zT

instance (ConvertibleExactly t1 t2) => (ConvertibleExactly (Complex t1) (Complex t2))
  where
  safeConvertExactly (a1 :+ i1) =
    do
    a2 <- safeConvertExactly a1
    i2 <- safeConvertExactly i1
    return $ a2 :+ i2

instance (HasEqAsymmetric a b) => HasEqAsymmetric (Complex a) (Complex b) where
  type EqCompareType (Complex a) (Complex b) = EqCompareType a b
  equalTo (a1 :+ i1) (a2 :+ i2) = (a1 == a2) && (i1 == i2)

instance (CanTestInteger t, CanTestZero t) => CanTestInteger (Complex t) where
  certainlyNotInteger (a :+ i) =
    certainlyNotInteger a || isCertainlyNonZero i
  certainlyIntegerGetIt (a :+ i) =
    case (certainlyIntegerGetIt a, certainlyIntegerGetIt i) of
      (Just aN, Just iN) | iN == 0 -> Just aN
      _ -> Nothing

instance CanNeg t => CanNeg (Complex t) where
  type NegType (Complex t) = Complex (NegType t)
  negate (a :+ i) = (negate a) :+ (negate i)

instance (CanAddAsymmetric a b) => CanAddAsymmetric (Complex a) (Complex b) where
  type AddType (Complex a) (Complex b) = Complex (AddType a b)
  add (a1 :+ i1) (a2 :+ i2) = (a1 + a2) :+ (i1 + i2)

instance (CanSub a b) => CanSub (Complex a) (Complex b) where
  type SubType (Complex a) (Complex b) = Complex (SubType a b)
  sub (a1 :+ i1) (a2 :+ i2) = (a1 - a2) :+ (i1 - i2)

instance
  (CanMulAsymmetric a b
  , CanAddSameType (MulType a b), CanSubSameType (MulType a b))
  =>
  CanMulAsymmetric (Complex a) (Complex b)
  where
  type MulType (Complex a) (Complex b) = Complex (MulType a b)
  mul (a1 :+ i1) (a2 :+ i2) =
    (a1*a2 - i1*i2) :+ (a1*i2 + i1*a2)

instance
  (CanMulAsymmetric a b
  , CanAddSameType (MulType a b), CanSubSameType (MulType a b)
  , CanMulAsymmetric b b, CanAddSameType (MulType b b)
  , CanDiv (MulType a b) (MulType b b))
  =>
  CanDiv (Complex a) (Complex b)
  where
  type DivType (Complex a) (Complex b) = Complex (DivType (MulType a b) (MulType b b))
  divide (a1 :+ i1) (a2 :+ i2) =
    let d = a2*a2 + i2*i2 in
    ((a1*a2 + i1*i2)/d) :+ ((i1*a2-a1*i2)/d)

instance
  (CanMulAsymmetric t t
  , CanAddSameType (MulType t t)
  , CanSqrt (MulType t t))
  =>
  CanAbs (Complex t)
  where
  type AbsType (Complex t) = SqrtType (MulType t t)
  abs (a :+ i) = sqrt (a*a + i*i)

instance
  (CanExp t
  , CanSinCos t
  , CanMulAsymmetric (ExpType t) (SinCosType t))
  =>
  CanExp (Complex t)
  where
  type ExpType (Complex t) = Complex (MulType (ExpType t) (SinCosType t))
  exp (a :+ i) =
    let ea = exp a in
    (ea * cos i) :+ (ea * sin i)

$(declForTypes
  [[t| Integer |], [t| Int |], [t| Rational |], [t| Double |]]
  (\ t -> [d|

    instance (HasEqAsymmetric $t b) => HasEqAsymmetric $t (Complex b) where
      type EqCompareType $t (Complex b) = EqCompareType $t b
      equalTo n (a2 :+ i2) = (n == a2) && (convertExactlyTargetSample n 0 == i2)

    instance (HasEqAsymmetric a $t) => HasEqAsymmetric (Complex a) $t where
      type EqCompareType (Complex a) $t = EqCompareType a $t
      equalTo (a1 :+ i1) n = (a1 == n) && (i1 == convertExactlyTargetSample n 0)

    instance (CanAddAsymmetric $t b) => CanAddAsymmetric $t (Complex b) where
      type AddType $t (Complex b) = Complex (AddType $t b)
      add n (a2 :+ i2) = (n + a2) :+ (convertExactlyTargetSample n 0 + i2)

    instance (CanAddAsymmetric a $t) => CanAddAsymmetric (Complex a) $t where
      type AddType (Complex a) $t = Complex (AddType a $t)
      add (a1 :+ i1) n = (a1 + n) :+ (i1 + (convertExactlyTargetSample n 0))

    instance (CanSub $t b) => CanSub $t (Complex b) where
      type SubType $t (Complex b) = Complex (SubType $t b)
      sub n (a2 :+ i2) = (n - a2) :+ (convertExactlyTargetSample n 0 - i2)

    instance (CanSub a $t) => CanSub (Complex a) $t where
      type SubType (Complex a) $t = Complex (SubType a $t)
      sub (a1 :+ i1) n = (a1 - n) :+ (i1 - (convertExactlyTargetSample n 0))

    instance
      (CanMulAsymmetric $t b) => CanMulAsymmetric $t (Complex b)
      where
      type MulType $t (Complex b) = Complex (MulType $t b)
      mul n (a2 :+ i2) = (n*a2) :+ (n*i2)

    instance
      (CanMulAsymmetric a $t) => CanMulAsymmetric (Complex a) $t
      where
      type MulType (Complex a) $t = Complex (MulType a $t)
      mul (a1 :+ i1) n = (a1*n) :+ (i1*n)

    instance
      (CanMulAsymmetric $t b
      , CanMulAsymmetric b b, CanAddSameType (MulType b b)
      , CanDiv (MulType $t b) (MulType b b))
      =>
      CanDiv $t (Complex b)
      where
      type DivType $t (Complex b) = Complex (DivType (MulType $t b) (MulType b b))
      divide n (a2 :+ i2) =
        let d = a2*a2 + i2*i2 in
        ((n*a2)/d) :+ (((-n)*i2)/d)

    instance
      (CanDiv a $t) => CanDiv (Complex a) $t
      where
      type DivType (Complex a) $t = Complex (DivType a $t)
      divide (a1 :+ i1) n = (a1/n) :+ (i1/n)
  |]))

max_stars_repo_name: michalkonecny/mixed-types-num

max_issues_repo_name: michalkonecny/mixed-types-num

max_forks_repo_name: michalkonecny/mixed-types-num


text:
{-# LANGUAGE ScopedTypeVariables #-}

module DecisionTree where

import qualified Data.Vector as V
import qualified Data.Vector.Generic as VG
import qualified Data.List as List
import qualified Data.List.Extras.Argmax as Argmax
import qualified Statistics.Sample as S

-- Returns the list of unique elements of x --
uniqueElems :: (Eq a) => V.Vector a -> [a]
uniqueElems x = List.nub $ V.toList x

-- Divide two Integral values --
intDiv :: (Integral a, Fractional b) => a -> a -> b
intDiv x y = fromIntegral x / fromIntegral y

-- Returns the sample probability mass function of `samples` P(X = x)
samplePMF :: (Eq a) => V.Vector a -> a -> Double
samplePMF samples x = V.length (V.filter (== x) samples) `intDiv` V.length samples

-- Returns the entropy from the samples `y`. In terms of the sample
-- frequencies p_j (for the j-the value) it is defined by
--                 H(Y) =  -\sum_j p_j * log_2 p_j
entropy :: forall a. (Eq a) => V.Vector a -> Double
entropy y = (-1.0) * sum (map partialEntropy (uniqueElems y))
    where partialEntropy :: a -> Double
          partialEntropy i = pr i * logBase 2 (pr i)

          pr :: a -> Double
          pr = samplePMF y

-- Calculates the "information gain" of the sample `y` when we subdivide
-- it into groups according to the labels given by `lab`
infoGain :: forall a b. (Eq a, Eq b) => V.Vector a -> V.Vector b -> Double
infoGain y lab
    | V.length y /= V.length lab = error $ "DecisionTree.hs:infoGain " ++
                                           "Dimensions do not match."
    | otherwise = (-1.0) * sum (map partialInfoGain subsets)
        where partialInfoGain :: V.Vector a -> Double
              partialInfoGain a = pr_a * entropy a
                  where pr_a = V.length a `intDiv` V.length y

              subsets :: [V.Vector a]
              subsets = map elemsWithLabel (uniqueElems lab)

              elemsWithLabel :: b -> V.Vector a
              elemsWithLabel a = VG.ifilter (\ i _ -> lab V.! i == a) y

fitTree :: forall a b. (Ord a, Eq b) => [V.Vector a] -> V.Vector b -> DNode [a] [Double]
fitTree x y = Decision ((decF bestP) . bestP)  (Result "left") (Result "right")
    where decF :: [a] -> a -> a -> Bool
          decF p z = z < S.mean (p x)

          projectors :: [[a] -> a]
          projectors = [\z -> (z !! n) | n <- [0..(length x - 1)]]
          bestP :: [a] -> a
          bestP = Argmax.argmax (\ p -> (infoGain y $ V.map (decF p) (p x))) projectors

-- a: experiment to decide on
-- b: type of result i.e. classification or classification probability
data DNode a b = Result b
               | Decision (a -> Bool) (DNode a b) (DNode a b)

-- Apply the decision tree to a single experiment
decide :: DNode a b -> a -> b
decide (Result p) = p
decide (Decision df tbranch fbranch) x =
        if df x
            then decide tbranch x
            else decide fbranch x

max_stars_repo_name: dseuss/stuff

max_issues_repo_name: dseuss/stuff

max_forks_repo_name: dseuss/stuff


text:
{-# LANGUAGE RecordWildCards #-}
module DFA where

import qualified Data.Set as Set
import qualified Data.Map.Strict as Map
import Control.Monad.State
import Data.Maybe
import Data.Semigroup (stimes)
import Numeric.LinearAlgebra hiding (remap, (<>))
import NFA
import Utils
import qualified Partition as P
import LruCache

-- | The DFA type is parameterized on the state type, so that we can use
-- Sets during the subset construction and then map them to integers for
-- better efficiency.
data DFA s c = DFA
  { dfaStart :: s
  , dfaFinal :: Set.Set s
  , dfaStates :: Set.Set s
  , dfaTransitions :: Set.Set (c, s, s)
  }
  deriving Show

-- | At first, the DFA states are sets of NFA states
type DfaState1 = Set.Set NfaState

-- | Then, the sets are mapped into integers
type DfaState2 = Int

----------------------------------------------------------------------
--                      The subset construction
----------------------------------------------------------------------

-- | Return all epsilon-edges from an NFA state
epsilonEdges :: Ord c => NFA (Maybe c) -> NfaState -> Set.Set NfaState
epsilonEdges NFA{nfaTransitions} s =
  fromMaybe mempty . Map.lookup Nothing . fromMaybe mempty . Map.lookup s $
    transitionMap nfaTransitions

-- | Compute an epsilon-closure
epsilonClosure :: Ord c => NFA (Maybe c) -> Set.Set NfaState -> Set.Set NfaState
epsilonClosure nfa = go mempty
  where
    go !seen0 !new0
      | Set.null new0 = seen0
      | otherwise =
      let
        seen1 = Set.union seen0 new0
        new1 = Set.difference (sconcatMap (epsilonEdges nfa) new0) seen0
      in go seen1 new1

-- | @move nfa s a@ is a set of NFA states where we can get from one of the
-- states in @s@ after receiving @a@ as an input
move :: Ord c => NFA (Maybe c) -> DfaState1 -> c -> DfaState1
move NFA{nfaTransitions} ss c = sconcatMap
  (\s -> fromMaybe mempty . Map.lookup (Just c) .
    fromMaybe mempty . Map.lookup s $
    transitionMap nfaTransitions)
  ss

-- | Build a DFA from an NFA
nfaToDfa :: forall c . Ord c => [c] -> NFA (Maybe c) -> DFA DfaState1 c
nfaToDfa alphabet nfa =
  let
    start :: DfaState1
    start = epsilonClosure nfa (nfaStart nfa)
    transitions = go Set.empty (Set.singleton $ start) Set.empty
    all_states = concat [ [s1,s2] | (_, s1, s2) <- Set.toList transitions ]
    contains_final :: DfaState1 -> Bool
    contains_final dfa_st = not . Set.null $ Set.intersection (nfaFinal nfa) dfa_st
  in
    DFA
      { dfaStart = start
      , dfaFinal = Set.fromList $ filter contains_final all_states
      , dfaStates = Set.fromList all_states
      , dfaTransitions = transitions
      }
  where
    go
      :: Set.Set DfaState1
      -> Set.Set DfaState1
      -> Set.Set (c, DfaState1, DfaState1)
      -> Set.Set (c, DfaState1, DfaState1)
    go !seen0 !new0 !transitions0
      | Set.null new0 = transitions0
      | otherwise =
      let
        seen1 = Set.union seen0 new0
        -- where can we get from the DFA states in new1?
        (transitions1, new1) =
          flip foldMap new0 $ \s ->
          flip foldMap alphabet $ \c ->
          let
            next :: DfaState1
            next = epsilonClosure nfa (move nfa s c)
          in (Set.singleton (c, s, next), Set.singleton next)

        transitions2 = transitions0 `Set.union` transitions1
        new2 = new1 `Set.difference` seen1
      in go seen1 new2 transitions2

----------------------------------------------------------------------
--                         DFA optimization
----------------------------------------------------------------------

-- | Map all DFA states
mapStates
  :: (Ord s1, Ord s2, Ord c)
  => (s1 -> s2)
  -> DFA s1 c
  -> DFA s2 c
mapStates remap DFA{..} =
  DFA
  { dfaStart = remap dfaStart
  , dfaFinal = Set.map remap dfaFinal
  , dfaStates = Set.map remap dfaStates
  , dfaTransitions = Set.map (\(c,s0,s1) -> (c, remap s0, remap s1)) dfaTransitions
  }

-- | Convert all states of the DFA to consecutive integers starting from 0
mapStatesToInt
  :: (Ord s, Ord c)
  => DFA s c
  -> DFA DfaState2 c
mapStatesToInt dfa =
  let
    stateMap = Map.fromList $ zip (Set.toList . dfaStates $ dfa) [0..]
    remap = (stateMap Map.!)
  in
    mapStates remap dfa

minimizeDfa
  :: forall s c . (Ord c, Ord s)
  => DFA s c
  -> DFA Int c
minimizeDfa dfa@DFA{..} = flip evalState 0 $ do
  initialPartition <- liftM2 (<>)
    (P.singleton dfaFinal)
    (P.singleton $ dfaStates `Set.difference` dfaFinal)

  let
    transitionMap :: Map.Map s (Map.Map c s)
    transitionMap = Map.fromListWith Map.union
      [ (s0, Map.singleton c s1)
      | (c, s0, s1) <- Set.toList dfaTransitions
      ]

    classify :: P.Partition s -> s -> Map.Map c Int
    classify pt s = fmap (P.lookupState pt) $ transitionMap Map.! s

  finalPartition <- P.subpartitionFixpoint classify initialPartition

  return $ mapStates (P.lookupState finalPartition) dfa

----------------------------------------------------------------------
--                Transfer matrices and probabilities
----------------------------------------------------------------------

data TransferMatrix = TransferMatrix
  { tmMatrix :: Matrix Double
    -- ^ transfer matrix itself
  , tmStart :: Vector Double
    -- ^ the indicator vector of the start state
  , tmFinal :: Vector Double
    -- ^ the indicator vector of the final states
  }
  deriving (Show)

-- | Construct the DFA transfer matrix assuming the uniform distribution
-- over nucleotides.
--
-- Return the matrix and the 0-based indices of the start and end states.
dfaTransferMatrix
  :: (Ord c, Ord s)
  => Map.Map c Double -- ^ character probabilities
  -> DFA s c
  -> TransferMatrix
dfaTransferMatrix freqs (mapStatesToInt -> DFA{..}) =
  -- we remap the DFA once again to make sure the numbers are consecutive
  let
    n = Set.size dfaStates
    elts =
      [ ((s0, s1), freqs Map.! c)
      | (c, s0, s1) <- Set.toList dfaTransitions
      ]
  in
    TransferMatrix
    { tmMatrix = (assoc (n,n) 0 . assocSum) elts
    , tmStart = assoc n 0 [ (dfaStart, 1) ]
    , tmFinal = assoc n 0 [ (st, 1) | st <- Set.toList dfaFinal ]
    }
  where
    assocSum = Map.toList . Map.fromListWith (+)

dfaProbability
  :: TransferMatrix
  -> Int -- ^ length of the random string where they motif may occur
  -> Double
dfaProbability TransferMatrix{..} len =
  (tmStart <# stimes len tmMatrix) <.> tmFinal

dfaProbabilityCached
  :: MonadLru Int (Matrix Double) m
  => TransferMatrix
  -> Int -- ^ length of the random string where they motif may occur
  -> m Double
dfaProbabilityCached TransferMatrix{..} len = do
  pw <- lruPow len tmMatrix
  return $ (tmStart <# pw) <.> tmFinal

max_stars_repo_name: feuerbach/motif-stats

max_issues_repo_name: feuerbach/motif-stats

max_forks_repo_name: feuerbach/motif-stats


text:
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
module Marvin.Test.TestUtils (
  isAround
  , nestedFromList
  , trainMatrix
  , equals
  , (+-)
) where

import Marvin.Test.Metric
import Marvin.API as Marvin

import Test.QuickCheck hiding (vector)
import qualified Test.QuickCheck as QC (vector)
import Numeric.LinearAlgebra.Data as LA
import Foreign.Storable

import qualified Data.Vector as Vec

nestedFromList = Vec.fromList . Prelude.map Vec.fromList

genVectorOfSize :: (Storable a, Arbitrary a) => Int -> Gen (Vector a)
genVectorOfSize n = do
  xs <- vectorOf n arbitrary
  return $ LA.fromList xs

instance (Storable a, Arbitrary a) => Arbitrary (Vector a) where
  arbitrary = sized genVectorOfSize

instance (Arbitrary a) => Arbitrary (Vec.Vector a) where
  arbitrary = Vec.fromList <$> arbitrary

instance Arbitrary NumericColumn where
  arbitrary = do
      t <- arbitrary
      h <- arbitrary
      let Right col = Marvin.fromList (h:t)
      return col

instance Arbitrary Error where
  arbitrary = do
    msg <- arbitrary
    return $ Marvin.Failure msg

isAround :: Metric a => a -> a -> Bool
isAround a x = dist a x <= threshold

infix 7 `equals`
equals :: Metric a => a -> a -> FloatPrecision -> Bool
equals x y (P e) = dist x y < e

newtype FloatPrecision = P Double

infix 8 +-
(+-) :: Double -> FloatPrecision
(+-) = P

threshold = 0.00001

trainMatrix = [
  [2,6,2,3,2,2,1,3,2,2,2,15.0],
  [1,6,1,3,3,1,1,2,2,4,2,14.5],
  [2,6,3,1,4,2,3,3,1,3,1,14.0],
  [2,6,1,4,4,2,2,3,2,4,2,16.0],
  [1,1,2,2,2,1,1,4,1,3,1,14.0],
  [1,3,2,1,1,1,2,2,1,2,1,16.5],
  [1,4,1,1,3,1,3,1,1,3,1,14.0],
  [1,3,2,3,2,1,2,3,1,3,1,16.5],
  [2,4,1,2,1,2,3,2,1,4,1,15.0],
  [2,4,1,1,2,2,3,3,1,3,1,16.0],
  [2,2,1,2,3,2,1,2,1,3,1,16.5],
  [2,5,2,2,1,2,3,1,2,1,2,15.0],
  [2,5,2,2,3,2,3,1,2,1,2,15.5],
  [2,6,1,4,3,2,2,3,2,4,2,15.5],
  [2,3,2,3,1,2,1,2,1,3,1,15.0],
  [1,1,2,4,3,1,1,1,1,2,1,15.0],
  [1,1,3,2,2,1,2,2,2,2,2,14.5],
  [1,2,1,1,2,1,2,1,1,4,1,16.0],
  [2,1,2,4,4,2,2,2,1,1,1,15.5],
  [2,2,2,2,3,2,3,4,1,2,1,14.5],
  [1,5,2,3,3,1,3,2,1,1,1,13.5],
  [2,5,3,4,4,2,2,4,2,3,2,16.0],
  [1,2,2,3,4,1,3,3,1,2,1,15.0],
  [2,3,1,2,2,2,2,2,1,1,1,14.5],
  [2,2,1,3,1,2,1,1,2,4,2,15.0],
  [1,1,2,3,4,1,1,2,1,4,1,14.5],
  [2,1,3,4,2,2,1,4,1,1,1,14.5],
  [1,3,3,4,3,1,3,4,2,4,2,15.5],
  [2,2,3,1,3,2,2,1,2,2,2,15.5],
  [1,5,1,2,4,1,2,1,2,4,2,15.5],
  [2,2,2,3,1,2,3,2,2,4,2,14.0],
  [1,2,3,4,2,1,1,4,1,3,1,17.0],
  [1,1,3,4,2,1,2,3,2,3,2,15.0],
  [1,2,1,2,1,1,2,3,1,3,1,16.0],
  [2,3,2,3,3,2,1,2,2,3,2,15.0],
  [2,5,1,1,4,2,1,4,2,1,2,16.0],
  [1,4,3,3,3,1,2,3,1,3,1,15.5],
  [1,2,3,3,3,1,1,3,2,1,2,16.5],
  [2,6,1,3,1,2,2,2,1,1,1,15.5],
  [1,6,2,1,1,1,2,3,2,3,2,15.0],
  [2,6,1,2,3,2,2,1,1,2,1,15.0],
  [2,2,1,4,2,2,1,4,2,1,2,15.0],
  [2,5,1,4,4,2,1,3,2,2,2,14.5],
  [1,5,1,3,1,1,2,2,2,2,2,14.5],
  [2,4,1,4,1,2,3,4,1,2,1,16.0],
  [1,6,1,4,2,1,1,4,2,3,2,14.5],
  [1,6,1,2,4,1,1,1,2,2,2,13.5],
  [2,6,1,1,2,2,2,4,1,3,1,15.0],
  [1,4,1,1,1,1,3,1,2,3,2,14.5],
  [2,1,1,4,3,2,3,3,2,3,2,15.0],
  [1,3,2,1,2,1,2,2,2,2,2,14.5],
  [2,5,1,2,2,2,1,1,1,4,1,15.5],
  [1,3,3,1,3,1,3,2,1,1,1,15.0],
  [1,2,3,4,3,1,1,4,2,3,2,16.5],
  [1,1,2,3,1,1,1,2,1,4,1,15.0],
  [1,5,3,1,2,1,1,1,1,2,1,14.0],
  [1,6,3,2,3,1,3,2,1,1,1,15.0],
  [2,5,3,1,3,2,2,3,1,4,1,14.0],
  [1,2,2,2,4,1,3,1,1,3,1,16.5],
  [2,1,2,3,4,2,2,3,1,4,1,16.0],
  [2,6,2,3,3,2,1,3,1,2,1,14.0],
  [1,3,3,1,1,1,3,2,2,1,2,15.5],
  [2,3,2,4,2,2,1,3,2,2,2,14.0],
  [1,5,2,2,1,1,3,4,1,4,1,14.5],
  [2,6,2,2,1,2,1,1,1,4,1,15.0],
  [1,2,3,1,4,1,1,2,2,4,2,16.0],
  [1,4,3,1,3,1,2,2,1,2,1,14.5],
  [2,5,3,4,3,2,2,4,2,3,2,15.0],
  [1,6,2,3,1,1,2,2,2,2,2,13.5],
  [2,4,1,3,4,2,3,1,1,1,1,15.0],
  [2,4,3,4,1,2,1,2,2,3,2,14.0],
  [1,4,3,4,1,1,2,4,2,1,2,15.5],
  [1,3,2,2,3,1,2,1,1,4,1,16.0],
  [1,6,3,3,2,1,3,4,1,4,1,15.0],
  [1,3,3,2,2,1,3,1,1,3,1,15.5],
  [1,6,2,4,3,1,2,1,2,4,2,15.5],
  [2,2,3,1,1,2,2,1,2,2,2,14.5],
  [2,2,2,3,2,2,3,2,1,4,1,14.5],
  [2,2,1,2,2,2,1,2,1,3,1,15.0],
  [2,5,2,3,2,2,3,3,1,3,1,15.0]]

testMatrix = [
  [1,1,1,4,4,1,3,4,1,4,1,16.0],
  [1,4,3,3,1,1,2,3,1,3,1,15.5],
  [2,2,1,4,4,2,1,4,2,1,2,15.0],
  [1,5,2,4,3,1,3,1,1,3,1,15.0],
  [1,4,2,1,2,1,1,2,2,4,2,14.5],
  [1,1,2,1,4,1,1,3,2,1,2,15.0],
  [1,1,3,1,3,1,2,1,2,4,2,15.0],
  [1,4,2,1,1,1,1,2,1,4,1,14.5],
  [2,6,1,4,2,2,2,3,1,4,1,15.0],
  [2,5,1,3,2,2,1,2,2,3,2,15.5],
  [1,6,1,3,1,1,1,2,1,4,1,15.0],
  [2,3,2,1,3,2,1,1,1,4,1,15.5],
  [1,2,2,3,3,1,3,3,2,2,2,14.0],
  [1,4,1,1,2,1,3,1,2,3,2,14.0],
  [1,2,1,2,4,1,2,3,2,3,2,15.0],
  [1,2,2,2,1,1,3,1,2,3,2,16.0],
  [1,5,1,1,1,1,2,3,1,3,1,15.5],
  [1,3,2,4,3,1,2,4,2,1,2,17.0],
  [2,1,3,4,3,2,1,4,2,1,2,15.5],
  [2,5,2,3,1,2,3,3,2,3,2,15.0]]

max_stars_repo_name: gaborhermann/marvin

max_issues_repo_name: gaborhermann/marvin

max_forks_repo_name: gaborhermann/marvin


text:
module Numerical.HBLAS.Lapack.FFI where
import Foreign.Ptr
import Foreign()
import Foreign.C.Types
import Data.Complex
import Data.Int

{-

stylenote: we will not use the LAPACKE_* operations, only the
LAPACKE_*_work variants that require an explicitly provided work buffer.

This is to ensure that solver routine allocation behavior is transparent

-}

{-
void LAPACK_dgesvx( char* fact, char* trans, lapack_int* n, lapack_int* nrhs,
                    double* a, lapack_int* lda, double* af, lapack_int* ldaf,
                    lapack_int* ipiv, char* equed, double* r, double* c,
                    double* b, lapack_int* ldb, double* x, lapack_int* ldx,
                    double* rcond, double* ferr, double* berr, double* work,
                    lapack_int* iwork, lapack_int *info );

-}

{-
    fortran FFI conventions!
-}

--type Stride_C =

newtype Fact_C = Fact_C CChar
newtype Trans_C = Trans_C CChar
newtype Stride_C = Stride_C Int32
newtype Equilib_C = Equilib_C CChar

type Fun_FFI_GESVX el = Ptr Fact_C  {- fact -}-> Ptr Trans_C {- trans -}
    -> Ptr Int32  {-n -}-> Ptr Int32 {- NRHS -}->
    Ptr el {- a -} -> Ptr Stride_C {- lda -} -> Ptr Double {- af -} -> Ptr Stride_C  {- ldf-}->
    Ptr Int32 -> Ptr Equilib_C {- equed -} -> Ptr el {- r -} -> Ptr el  ->
    Ptr el {- b -} -> Ptr Stride_C {- ld b   -} -> Ptr el {- x -} -> Ptr Stride_C {- ldx -}->
    Ptr el {-rcond -}-> Ptr el {- ferr-} -> Ptr el {-berr-} -> Ptr el {-work-}->
    Ptr Int32 {-iwork -}-> Ptr Int32 {-info  -} -> IO ()

{-

the prefixes mean s=single,d=double,c=complex float,d=complex double

fact will be a 1 character C string
either
    "F", then the inputs af and ipiv already contain the permuted LU factorization
        (act as input rather than result params)
    "E", Matrix input A will be equilibriated if needed, then copied to AF and Factored
    "N", matrix input A will be copied to AF

-}

{-
Xgesvx  is the s -sing

im assuming for now that any real use of *gesvx routines, or any other
n^3 complexity algs from LAPACK, are on inputs typically  n>=15, which means > 1000 flops,
which is > 1µs, and thus ok to
-}

--need to get around to wrapping these, but thats for another day
foreign import ccall  "sgesvx_"  sgesvx :: Fun_FFI_GESVX Float
foreign import ccall  "dgesvx_"  dgesvx :: Fun_FFI_GESVX Double
foreign import ccall  "cgesvx_"  cgesvx :: Fun_FFI_GESVX (Complex Float)
foreign import ccall  "zgesvx_"  zgesvx :: Fun_FFI_GESVX (Complex Double)

--lapack_int ?syev_(  char *jobz, char *uplo, lapack_int *n, ?* a, lapack_int * lda, ?* w );
-- ? is Double or Float

newtype JobTy = JBT CChar
newtype UploTy = UPLT CChar
newtype Info = Info Int32

--basic symmetric eigen value solvers
type SYEV_FUN_FFI elem = Ptr JobTy -> Ptr UploTy -> Ptr Int32  -> Ptr elem -> Ptr Int32 -> Ptr elem -> Ptr Info-> IO ()
foreign import ccall "ssyev_" ssyev_ffi :: SYEV_FUN_FFI Float
foreign import ccall "dsyev_" dsyev_ffi :: SYEV_FUN_FFI Double

{-unsafe versions of lapack routines are meant to ONLY be used for workspace queries-}
foreign import ccall unsafe "ssyev_" ssyev_ffi_unsafe :: SYEV_FUN_FFI Float
foreign import ccall unsafe "dsyev_" dsyev_ffi_unsafe :: SYEV_FUN_FFI Double

--lapack_int LAPACKE_<?>gesv( int matrix_order, lapack_int n, lapack_int nrhs, <datatype>* a, lapack_int lda, lapack_int* ipiv, <datatype>* b, lapack_int ldb );
--call sgesv( n, nrhs, a, lda, ipiv, b, ldb, info )

type GESV_FUN_FFI elem = Ptr Int32  {- n -} -> Ptr Int32 {- nrhs -} -> Ptr elem {- a -}
            -> Ptr Int32 {-  lda -} -> Ptr Int32 {- permutation vector -}
            -> Ptr elem  {- b -} -> Ptr Int32 {- ldb -} -> Ptr Info -> IO ()
-- | basic Linear system solvers. they act inplace
foreign import ccall  "sgesv_"  sgesv_ffi  ::GESV_FUN_FFI Float
foreign import ccall  "dgesv_"  dgesv_ffi :: GESV_FUN_FFI Double
foreign import ccall  "cgesv_"  cgesv_ffi :: GESV_FUN_FFI (Complex Float)
foreign import ccall  "zgesv_"  zgesv_ffi :: GESV_FUN_FFI (Complex Double)

-- / not sure if linear solvers ever are run in < 1 microsecond size instances in practice
foreign import ccall unsafe "sgesv_"  sgesv_ffi_unsafe  ::GESV_FUN_FFI Float
foreign import ccall unsafe "dgesv_"  dgesv_ffi_unsafe :: GESV_FUN_FFI Double
foreign import ccall unsafe "cgesv_"  cgesv_ffi_unsafe :: GESV_FUN_FFI (Complex Float)
foreign import ccall unsafe "zgesv_"  zgesv_ffi_unsafe :: GESV_FUN_FFI (Complex Double)

{- only provide  -}

max_stars_repo_name: archblob/hblas

max_issues_repo_name: archblob/hblas

max_forks_repo_name: archblob/hblas


text:
module Mnist.Load where

  import Neural
  import Numeric.LinearAlgebra
  import Debug.Trace
  import qualified Data.ByteString.Lazy as B
  import Data.Binary.Get
  import Data.Int
  import Control.Monad

  -- MNIST binary data parsers
  getLabelHeader :: Get Int
  getLabelHeader = do
    _ <- getWord32be
    num <- getInt32be
    return $ fromIntegral num

  getLabel :: Get (Vector Double)
  getLabel = do
    label <- getInt8
    return $ (toOneHot . fromIntegral) label
      where
        toOneHot :: Int -> Vector Double
        toOneHot v = vector $ (\i -> if i==v then 1.0 else 0.0) `map` [0..9]

  genVector :: Get (Vector Double) -> Get [Vector Double]
  genVector getOne = do
    empty <- isEmpty
    if empty then return []
    else do h <- getOne
            tl <- genVector getOne
            return (h:tl)

  getImageHeader :: Get (Int, Int, Int)
  getImageHeader = do
    _ <- getWord32be
    num <- getInt32be
    w <- getInt32be
    h <- getInt32be
    return (fromIntegral num, fromIntegral w, fromIntegral h)

  getImage :: Get (Vector Double)
  getImage = do
    colors <- replicateM (28*28) getWord8
    return $ vector ((toRealColor . fromIntegral) `map` colors)

  toRealColor :: Int -> Double
  toRealColor v = (fromIntegral v) / 255

  -- Parse MNIST Dataset into a list of examples (pair of a list of images and a list of labels)
  loadExamples :: IO [Example]
  loadExamples = do
    rawImage1 <- B.readFile "data/train-images.idx3-ubyte"
    rawImage2 <- B.readFile "data/t10k-images.idx3-ubyte"
    rawLabel1 <- B.readFile "data/train-labels.idx1-ubyte"
    rawLabel2 <- B.readFile "data/t10k-labels.idx1-ubyte"
    let labels1 = runGet (getLabelHeader *> (genVector getLabel)) rawLabel1
    let labels2 = runGet (getLabelHeader *> (genVector getLabel)) rawLabel2
    let images1 = runGet (getImageHeader *> (genVector getImage)) rawImage1
    let images2 = runGet (getImageHeader *> (genVector getImage)) rawImage2
    return $ zip (images1++images2) (labels1++labels2)

max_stars_repo_name: ycjungSubhuman/hskmlp

max_issues_repo_name: ycjungSubhuman/hskmlp

max_forks_repo_name: ycjungSubhuman/hskmlp


text:
{-# LANGUAGE FlexibleContexts #-}

import Data.Function
import Data.Fixed
import Numeric.LinearAlgebra

grid :: Integer -> Matrix Double
grid n = build (300, 300) (builder `on` (+1) . round)
  where builder x y = let id = x + 10
                       in fromInteger . subtract 5
                        $ (id * y + n) * id `mod` 1000 `div` 100

main = getContents >>= putStrLn
                     . (\(x, y) -> show (x+1) ++ "," ++ show (y+1))
                     . maxIndex
                     . corr2 (konst 1 (3, 3))
                     . grid . read