text:
{-# LANGUAGE ViewPatterns #-}
-- | Tests for Statistics.Math
module Tests.SpecFunctions (
  tests
  ) where

import qualified Data.Vector as V
import           Data.Vector   ((!))

import Test.QuickCheck  hiding (choose,within)
import Test.Framework
import Test.Framework.Providers.QuickCheck2
import Test.Framework.Providers.HUnit
import Test.HUnit (assertBool)

import Tests.Helpers
import Tests.SpecFunctions.Tables
import Numeric.SpecFunctions
import Numeric.MathFunctions.Comparison (within,relativeError)
import Numeric.MathFunctions.Constants  (m_epsilon,m_tiny)

tests :: Test
tests = testGroup "Special functions"
  [ testProperty "Gamma(x+1) = x*Gamma(x) [logGamma]"  $ gammaReccurence logGamma  3e-8
  , testProperty "Gamma(x+1) = x*Gamma(x) [logGammaL]" $ gammaReccurence logGammaL 2e-13
  , testProperty "gamma(1,x) = 1 - exp(-x)"      $ incompleteGammaAt1Check
  , testProperty "0 <= gamma <= 1"               $ incompleteGammaInRange
  , testProperty "0 <= I[B] <= 1"            $ incompleteBetaInRange
  , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse
  -- XXX FIXME DISABLED due to failures
  -- , testProperty "invIncompleteBeta  = B^-1" $ invIBetaIsInverse
  , testProperty "gamma - increases" $
      \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y
  , testProperty "invErfc = erfc^-1"         $ invErfcIsInverse
  , testProperty "invErf  = erf^-1"          $ invErfIsInverse
    -- Unit tests
  , testAssertion "Factorial is expected to be precise at 1e-15 level"
      $ and [ eq 1e-15 (factorial (fromIntegral n :: Int))
                       (fromIntegral (factorial' n))
            |n <- [0..170]]
  , testAssertion "Log factorial is expected to be precise at 1e-15 level"
      $ and [ eq 1e-15 (logFactorial (fromIntegral n :: Int))
                       (log $ fromIntegral $ factorial' n)
            | n <- [2..170]]
  , testAssertion "logGamma is expected to be precise at 1e-9 level [integer points]"
      $ and [ eq 1e-9 (logGamma (fromIntegral n))
                      (logFactorial (n-1))
            | n <- [3..10000::Int]]
  , testAssertion "logGamma is expected to be precise at 1e-9 level [fractional points]"
      $ and [ eq 1e-9 (logGamma x) lg | (x,lg) <- tableLogGamma ]
  , testAssertion "logGammaL is expected to be precise at 1e-15 level"
      $ and [ eq 1e-15 (logGammaL (fromIntegral n))
                       (logFactorial (n-1))
            | n <- [3..10000::Int]]
  , testAssertion "logGammaL is expected to be precise at 1e-10 level [fractional points]"
      $ and [ eq (64*m_epsilon) (logGammaL x) lg | (x,lg) <- tableLogGamma ]
    -- FIXME: loss of precision when logBeta p q ≈ 0.
    --        Relative error doesn't work properly in this case.
  , testAssertion "logBeta is expected to be precise at 1e-6 level"
      $ and [ eq 1e-6 (logBeta p q)
                      (logGammaL p + logGammaL q - logGammaL (p+q))
            | p <- [0.1,0.2 .. 0.9] ++ [2 .. 20]
            , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]
            ]
  , testAssertion "digamma is expected to be precise at 1e-14 [integers]"
      $ digammaTestIntegers 1e-14
    -- Relative precision is lost when digamma(x) ≈ 0
  , testAssertion "digamma is expected to be precise at 1e-12"
      $ and [ eq 1e-12 r (digamma x) | (x,r) <- tableDigamma ]
    --
  , let deviations = [ ( "p=",p, "q=",q, "x=",x
                       , "ib=",ib, "ib'=",ib'
                       , "err=",relativeError ib ib' / m_epsilon)
                     | (p,q,x,ib) <- tableIncompleteBeta
                     , let ib' = incompleteBeta p q x
                     , not $ eq (64 * m_epsilon) ib' ib
                     ]
    in testCase "incompleteBeta is expected to be precise at 32*m_epsilon level"
     $ assertBool (unlines (map show deviations)) (null deviations)
  , testAssertion "incompleteBeta with p > 3000 and q > 3000"
      $ and [ eq 1e-11 (incompleteBeta p q x) ib | (x,p,q,ib) <-
                 [ (0.495,  3001,  3001, 0.2192546757957825068677527085659175689142653854877723)
                 , (0.501,  3001,  3001, 0.5615652382981522803424365187631195161665429270531389)
                 , (0.531,  3500,  3200, 0.9209758089734407825580172472327758548870610822321278)
                 , (0.501, 13500, 13200, 0.0656209987264794057358373443387716674955276089622780)
                 ]
            ]
  , testAssertion "choose is expected to precise at 1e-12 level"
      $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k)
            | n <- [0..1000], k <- [0..n]]
  , testAssertion "logChoose == log . choose"
      $ and [ let n' = fromIntegral n
                  k' = fromIntegral k
              in within 2 (logChoose n' k') (log $ choose n' k')
            | n <- [0::Int .. 1000], k <- [0 .. n]]
    ----------------------------------------------------------------
    -- Self tests
  , testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0
  ]

----------------------------------------------------------------
-- QC tests
----------------------------------------------------------------

-- Γ(x+1) = x·Γ(x)
gammaReccurence :: (Double -> Double) -> Double -> Double -> Property
gammaReccurence logG ε x =
  (x > 0 && x < 100)  ==>  (abs (g2 - g1 - log x) < ε)
    where
      g1 = logG x
      g2 = logG (x+1)

-- γ(s,x) is in [0,1] range
incompleteGammaInRange :: Double -> Double -> Property
incompleteGammaInRange (abs -> s) (abs -> x) =
  x >= 0 && s > 0  ==> let i = incompleteGamma s x in i >= 0 && i <= 1

-- γ(1,x) = 1 - exp(-x)
-- Since Γ(1) = 1 normalization doesn't make any difference
incompleteGammaAt1Check :: Double -> Property
incompleteGammaAt1Check (abs -> x) =
  x > 0 ==> (incompleteGamma 1 x + exp(-x)) ≈ 1
  where
    (≈) = eq 1e-13

-- invIncompleteGamma is inverse of incompleteGamma
invIGammaIsInverse :: Double -> Double -> Property
invIGammaIsInverse (abs -> a) (range01 -> p) =
  a > m_tiny && p > m_tiny && p < 1  ==>
    ( counterexample ("a    = " ++ show a )
    $ counterexample ("p    = " ++ show p )
    $ counterexample ("x    = " ++ show x )
    $ counterexample ("p'   = " ++ show p')
    $ counterexample ("err  = " ++ show (relativeError p p'))
    $ counterexample ("pred = " ++ show δ)
    $ relativeError p p' < δ
    )
  where
    x  = invIncompleteGamma a p
    f' = exp ( log x * (a-1) - x - logGamma a)
    p' = incompleteGamma    a x
    -- FIXME: 128 is big constant. It should be replaced by something
    --        smaller when #42 is fixed
    δ  = (m_epsilon/2) * (256 + 1 * (1 + abs (x * f' / p)))

-- invErfc is inverse of erfc
invErfcIsInverse :: Double -> Property
invErfcIsInverse ((*2) . range01 -> p)
  = counterexample ("p  = " ++ show p )
  $ counterexample ("x  = " ++ show x )
  $ counterexample ("p' = " ++ show p')
  $ abs (p - p') <= 1e-14
  where
    x  = invErfc p
    p' = erfc x

-- invErf is inverse of erf
invErfIsInverse :: Double -> Property
invErfIsInverse a
  = counterexample ("p  = " ++ show p )
  $ counterexample ("x  = " ++ show x )
  $ counterexample ("p' = " ++ show p')
  $ abs (p - p') <= 1e-14
  where
    x  = invErf p
    p' = erf x
    p  | a < 0     = - range01 a
       | otherwise =   range01 a

-- B(s,x) is in [0,1] range
incompleteBetaInRange :: Double -> Double -> Double -> Property
incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =
  p > 0 && q > 0  ==> let i = incompleteBeta p q x in i >= 0 && i <= 1

-- invIncompleteBeta is inverse of incompleteBeta
invIBetaIsInverse :: Double -> Double -> Double -> Property
invIBetaIsInverse (abs -> p) (abs -> q) (range01 -> x) =
  p > 0 && q > 0  ==> ( counterexample ("p   = " ++ show p )
                      $ counterexample ("q   = " ++ show q )
                      $ counterexample ("x   = " ++ show x )
                      $ counterexample ("x'  = " ++ show x')
                      $ counterexample ("a   = " ++ show a)
                      $ counterexample ("err = " ++ (show $ abs $ (x - x') / x))
                      $ abs (x - x') <= 1e-12
                      )
  where
    x' = incompleteBeta    p q a
    a  = invIncompleteBeta p q x

-- Table for digamma function:
--
-- Uses equality ψ(n) = H_{n-1} - γ where
--   H_{n} = Σ 1/k, k = [1 .. n]     - harmonic number
--   γ     = 0.57721566490153286060  - Euler-Mascheroni number
digammaTestIntegers :: Double -> Bool
digammaTestIntegers eps
  = all (uncurry $ eq eps) $ take 3000 digammaInt
  where
    ok approx exact = approx
    -- Harmonic numbers starting from 0
    harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]
    gam   = 0.57721566490153286060
    -- Digamma values
    digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN

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

-- Lookup table for fact factorial calculation. It has fixed size
-- which is bad but it's OK for this particular case
factorial_table :: V.Vector Integer
factorial_table = V.generate 2000 (\n -> product [1..fromIntegral n])

-- Exact implementation of factorial
factorial' :: Integer -> Integer
factorial' n = factorial_table ! fromIntegral n

-- Exact albeit slow implementation of choose
choose' :: Integer -> Integer -> Integer
choose' n k = factorial' n `div` (factorial' k * factorial' (n-k))

-- Truncate double to [0,1]
range01 :: Double -> Double
range01 = abs . (snd :: (Integer, Double) -> Double) . properFraction

max_stars_repo_name: Javran/math-functions

max_issues_repo_name: Javran/math-functions

max_forks_repo_name: Javran/math-functions


text:
module SchemeParser.Parser where

import Numeric (readInt, readFloat, readHex, readOct)
import Data.Ratio
import Data.Complex
import Data.Char (digitToInt)
import Data.Maybe (listToMaybe, fromJust)

import qualified Data.Map as M

import Data.Functor (($>))

import Text.Parsec hiding (spaces)
import Text.Parsec.Prim

import SchemeParser.Types

-- parser helpers

symbol :: LispParser Char
symbol = oneOf "!$%&|*+-/:<=>?@^_~"

spaces :: LispParser ()
spaces = skipMany1 space

escapedChars :: LispParser Char
escapedChars = do
  char '\\'
  c <- oneOf "\"nrt\\"
  return $ case c of
    'n' -> '\n'
    'r' -> '\r'
    't' -> '\t'
    _   -> c

readBin :: Integral a => String -> Maybe a
readBin = fmap fst . listToMaybe . readInt 2 (`elem` "01") digitToInt

--
-- Lisp value parsers
--

-- primitives

parseAtom :: LispParser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many (letter <|> digit <|> symbol)
  let atom = first : rest
  return $ LAtom atom

parseString :: LispParser LispVal
parseString = do
  char '"'
  x <- many $ escapedChars <|> noneOf "\""
  char '"'
  return $ LString x

parseChar :: LispParser LispVal
parseChar = do
  try (string "#\\")
  c <- try (string "newline" <|> string "space")
       <|> do { x <- anyChar; notFollowedBy alphaNum; return [x]}
  return $ LChar $ case c of
    "newline" -> '\n'
    "space"   -> ' '
    _         -> head c

parseBool :: LispParser LispVal
parseBool =
  try (string "#t" $> LBool True) <|> try (string "#f" $> LBool False)

-- http://www.schemers.org/Documents/Standards/R5RS/HTML/r5rs-Z-H-9.html#%_sec_6.2.4
parseNumberRadix :: LispParser LispVal
parseNumberRadix =
  parseBin <|> parseOct <|> parseDec <|> parseHex
  where
    parseBin = do
      try (string "#b")
      LNumber . fromJust . readBin <$> many1 (oneOf "01")
    parseOct = do
      try (string "#o")
      LNumber . fst . (!! 0) . readOct <$> many1 octDigit
    parseDec = do
      try (string "#d")
      parseNumber
    parseHex = do
      try (string "#x")
      LNumber . fst . (!! 0) . readHex <$> many1 hexDigit

parseNumber :: LispParser LispVal
parseNumber = parseNumberRadix <|> LNumber . read <$> many1 digit

parseFloat :: LispParser LispVal
parseFloat = do
  l <- many1 digit
  char '.'
  r <- many1 digit
  return $ LFloat $ fst $ head $ readFloat (l ++ "." ++ r)

parseRatio :: LispParser LispVal
parseRatio = do
  n <- many1 digit
  char '/'
  d <- many1 digit
  return $ LRatio (read n % read d)

parseComplex :: LispParser LispVal
parseComplex = do
  r <- try parseFloat <|> parseNumber
  char '+'
  i <- try parseFloat <|> parseNumber
  char 'i'
  return $ LComplex (toDouble r :+ toDouble i)

-- lists

parseList :: LispParser LispVal
parseList = LList <$> sepBy parseExpr spaces

parseDottedList :: LispParser LispVal
parseDottedList = do
  hd <- endBy parseExpr spaces
  tl <- char '.' >> spaces >> parseExpr
  return $ LDottedList hd tl

parseQuoted :: LispParser LispVal
parseQuoted = do
  char '\''
  x <- parseExpr
  return $ LList [LAtom "quote", x]

parseQuasiQuoted :: LispParser LispVal
parseQuasiQuoted = do
  char '`'
  x <- parseExpr
  return $ LList [LAtom "quasiquote", x]

parseUnquoted :: LispParser LispVal
parseUnquoted = do
  char ','
  x <- parseExpr
  return $ LList [LAtom "unquote", x]

parseUnquoteSpliced :: LispParser LispVal
parseUnquoteSpliced = do
  try $ string ",@"
  x <- parseExpr
  return $ LList [LAtom "unquote-splice", x]

parseVector :: LispParser LispVal
parseVector = LVector <$> sepBy parseExpr spaces

--
-- helpers
--

toDouble :: LispVal -> Double
toDouble (LFloat f) = realToFrac f
toDouble (LNumber n) = fromIntegral n

--

parseExpr :: LispParser LispVal
parseExpr = parseAtom
            <|> parseString
            <|> parseChar
            <|> try parseFloat
            <|> try parseRatio
            <|> try parseComplex
            <|> parseNumber
            <|> parseBool
            <|> parseQuoted
            <|> parseQuasiQuoted
            <|> parseUnquoteSpliced
            <|> parseUnquoted
            <|> char '(' *> (try parseList <|> parseDottedList) <* char ')'
            <|> try (string "#(" *> parseVector <* char ')')

ext: hs

max_stars_repo_name: alexpeits/haskell-lab-old

max_issues_repo_name: alexpeits/haskell-lab-old

max_forks_repo_name: alexpeits/haskell-lab-old


text:
module Naive
(
  reduceByRowsV
  , reduceByColumnsV
)

where

import qualified Numeric.LinearAlgebra as LA

type Matrix = LA.Matrix LA.R
type Vector = LA.Vector LA.R
type R = LA.R

reduceByRowsV :: (Vector -> R) -> Matrix -> Vector
reduceByRowsV f = LA.vector . map f . LA.toRows

reduceByColumnsV :: (Vector -> R) -> Matrix -> Vector
reduceByColumnsV f = LA.vector . map f . LA.toColumns

max_stars_repo_name: Alexander-Ignatyev/morpheus

max_issues_repo_name: aligusnet/morpheus

max_forks_repo_name: aligusnet/morpheus


text:
module NonlinearTransform
  ( solution8
  , solution9
  , solution10
  , createTransformedX2
  , createTransformedX
  , createTransformedXk
  ) where

import           LinearRegression
--------------------------------------------------------------
-- Solutions to homework 2 of "Learning from data" question 8-10
--------------------------------------------------------------
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.HMatrix

-- | the target function
targetf :: (R, R) -> R
targetf (x1, x2) = signum (x1 ^ 2 + x2 ^ 2 - 0.6)

-- | returns the second value in a 3-tuple
sndOfTriple :: (a, b, c) -> b
sndOfTriple (_, b, _) = b

-- | returns the first value in a 3-tuple
fstOfTriple :: (a, b, c) -> a
fstOfTriple (a, _, _) = a

-- | Creates the transformed matrix X with feature vectors (1,x1,x2,x1*x2,x1^2,x2^2) from a list of data points
createTransformedX :: Int -> ([(R, R)] -> [R]) -> [(R, R)] -> Matrix R
createTransformedX numFeatures transformfunction listOfPoints =
  matrix numFeatures listOfNumbers
  where
    listOfNumbers = transformfunction listOfPoints

transformFunction1 =
  \listOfPoints ->
    concat [[1, a, b, a * b, a ^ 2, b ^ 2] | (a, b) <- listOfPoints]

transformFunction2 =
  \listOfPoints ->
    concat
      [ [1, a, b, a * b, a ^ 2, b ^ 2, abs (a - b), abs (a + b)]
      | (a, b) <- listOfPoints
      ]

transformFunctionk k =
  \listOfPoints ->
    concat
      [ take (k + 1) [1, a, b, a ^ 2, b ^ 2, a * b, abs (a - b), abs (a + b)]
      | (a, b) <- listOfPoints
      ]

createTransformedX1 listOfPoints =
  createTransformedX 6 transformFunction1 listOfPoints

createTransformedX2 listOfPoints =
  createTransformedX 8 transformFunction2 listOfPoints

createTransformedXk k listOfPoints =
  createTransformedX (k + 1) (transformFunctionk k) listOfPoints

-- | solution to Problem 8
solution8 :: IO R
solution8 = fmap sndOfTriple $ trainLinRegWNoise targetf 1000 0.1 createMatrixX

-- | solution to Problem 9
solution9 :: IO (Vector R)
solution9 =
  fmap fstOfTriple $ trainLinRegWNoise targetf 1000 0.1 createTransformedX1

-- | solution to Problem 10
solution10 :: IO R
solution10 = solution9 >>= testLinRegWNoise targetf 1000 0.1 createTransformedX1


text:
{-# LANGUAGE DeriveFunctor #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE ImportQualifiedPost #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns #-}
{-# OPTIONS_GHC -Wwarn #-}

module Data.Distribution
  ( ToRealFrac(..)
  , Distribution(..)
  , computeDistribution
  , computeDistributionStats
  , mapToDistribution
  , zeroDistribution
  , dPercIx
  , dPercSpec
  , dPercSpec'
  , PercSpec(..)
  , renderPercSpec
  , Percentile(..)
  , pctFrac
  , stdPercentiles
  -- Aux
  , spans
  ) where

import Prelude (String, (!!), error, head, last, show)
import Cardano.Prelude hiding (head, show)

import Control.Arrow
import Data.Aeson (FromJSON(..), ToJSON(..))
import Data.Foldable qualified as F
import Data.List (span)
import Data.Vector (Vector)
import Data.Vector qualified as Vec
import Statistics.Sample qualified as Stat
import Text.Printf (PrintfArg, printf)

data Distribution a b =
  Distribution
  { dSize         :: Int
  , dAverage      :: a
  , dPercentiles  :: [Percentile a b]
  }
  deriving (Functor, Generic, Show)

instance (FromJSON a, FromJSON b) => FromJSON (Distribution a b)
instance (  ToJSON a,   ToJSON b) => ToJSON   (Distribution a b)

newtype PercSpec a = Perc { psFrac :: a } deriving (Eq, Generic, Show)

instance (FromJSON a) => FromJSON (PercSpec a)
instance (  ToJSON a) => ToJSON   (PercSpec a)

dPercIx :: Int -> Distribution a b -> b
dPercIx i d = pctSample $ dPercentiles d !! i

dPercSpec :: Eq (PercSpec a) => PercSpec a -> Distribution a b -> Maybe b
dPercSpec p = fmap pctSample . find ((== p) . pctSpec) . dPercentiles

dPercSpec' :: (Show a, Eq (PercSpec a)) => String -> PercSpec a -> Distribution a b -> b
dPercSpec' desc p =
  maybe (error er) pctSample . find ((== p) . pctSpec) . dPercentiles
 where
   er = printf "Missing centile %f in %s" (show $ psFrac p) desc

renderPercSpec :: PrintfArg a => Int -> PercSpec a -> String
renderPercSpec width = \case
  Perc x    -> printf ("%0."<>show (width-2)<>"f") x

data Percentile a b =
  Percentile
  { pctSpec        :: !(PercSpec a)
  , pctSample      :: !b
  }
  deriving (Functor, Generic, Show)

instance (FromJSON a, FromJSON b) => FromJSON (Percentile a b)
instance (  ToJSON a,   ToJSON b) => ToJSON   (Percentile a b)

pctFrac :: Percentile a b -> a
pctFrac = psFrac . pctSpec

stdPercentiles :: [PercSpec Float]
stdPercentiles =
  [ Perc 0.01, Perc 0.05
  , Perc 0.1, Perc 0.2, Perc 0.3, Perc 0.4
  , Perc 0.5, Perc 0.6
  , Perc 0.7, Perc 0.75
  , Perc 0.8, Perc 0.85, Perc 0.875
  , Perc 0.9, Perc 0.925, Perc 0.95, Perc 0.97, Perc 0.98, Perc 0.99
  , Perc 0.995, Perc 0.997, Perc 0.998, Perc 0.999
  , Perc 0.9995, Perc 0.9997, Perc 0.9998, Perc 0.9999
  ]

zeroDistribution :: Num a => Distribution a b
zeroDistribution =
  Distribution
  { dSize        = 0
  , dAverage     = 0
  , dPercentiles = mempty
  }

-- | For a list of distributions, compute a distribution of averages and rel stddev
-- (aka. coefficient of variance).
computeDistributionStats ::
    forall a v
  . ( v ~ Double -- 'v' is fixed by Stat.stdDev
    , Num a
    )
  => String -> [Distribution a v]
  -> Either String (Distribution a v, Distribution a v)
computeDistributionStats desc xs = do
  when (null xs) $
    Left $ "Empty list of distributions in " <> desc
  let distPcts    = dPercentiles <$> xs
      pctDistVals = transpose distPcts
  unless (all (pctLen ==) (length <$> distPcts)) $
    Left ("Distributions with different percentile counts: " <> show (length <$> distPcts) <> " in " <> desc)
  pure $ (join (***) (Distribution (length xs) 0)
          :: ([Percentile a v], [Percentile a v]) -> (Distribution a v, Distribution a v))
       $ unzip (pctsMeanCoV <$> pctDistVals)
 where
   pctLen = length . dPercentiles $ head xs

   pctsMeanCoV :: [Percentile a v] -> (Percentile a v, Percentile a v)
   pctsMeanCoV xs' = join (***) (Percentile . pctSpec $ head xs')
     (mean, Stat.stdDev vec / mean)
    where
      vec = Vec.fromList $ pctSample <$> xs'
      mean = Stat.mean vec

mapToDistribution :: (Real v, ToRealFrac v a) => (b -> v) -> [PercSpec a] -> [b] -> Distribution a v
mapToDistribution f pspecs xs = computeDistribution pspecs (f <$> xs)

computeDistribution :: (Real v, ToRealFrac v a) => [PercSpec a] -> [v] -> Distribution a v
computeDistribution percentiles (sort -> sorted) =
  Distribution
  { dSize        = size
  , dAverage     = toRealFrac (F.sum sorted) / fromIntegral (size `max` 1)
  , dPercentiles =
    (Percentile     (Perc 0)   mini:) .
    (<> [Percentile (Perc 1.0) maxi]) $
    percentiles <&>
      \spec ->
        let sample = if size == 0
                     then 0
                     else sorted !! indexAtFrac (psFrac spec)
        in Percentile spec sample
  }
 where size = length sorted
       indexAtFrac f = floor (fromIntegral (size - 1) * f)
       (,) mini maxi =
         if size == 0
         then (0,           0)
         else (head sorted, last sorted)

class RealFrac b => ToRealFrac a b where
  toRealFrac :: a -> b

instance RealFrac b => ToRealFrac Int b where
  toRealFrac = fromIntegral

instance {-# OVERLAPPABLE #-} (RealFrac b, Real a) => ToRealFrac a b where
  toRealFrac = realToFrac

spans :: forall a. (a -> Bool) -> [a] -> [Vector a]
spans f = go []
 where
   go :: [Vector a] -> [a] -> [Vector a]
   go acc [] = reverse acc
   go acc xs =
     case span f $ dropWhile (not . f) xs of
       ([], rest) -> go acc rest
       (ac, rest) ->
         go (Vec.fromList ac:acc) rest

max_stars_repo_name: Adlistic/cardano-node

max_issues_repo_name: Adlistic/cardano-node

max_forks_repo_name: Adlistic/cardano-node


text:
{-| Bindings to FFTW.

Example usage:

> import Foreign.Marshal.Array
> import Data.Complex
> import Foreign.Storable.Complex
> import FFTW
>
> main = do
>     inA  <- fftwAllocComplex 1024
>     outA <- fftwAllocComplex 1024
>
>     plan <- planDFT1d 1024 inA outA Forward fftwEstimate
>
>     pokeArray inA $ map (:+ 0) [0..1023]
>     execute plan
>     res <- peekArray 1024 outA
>
>     fftwFree inA
>     fftwFree outA
>
>     print res
-}

module Numeric.FFTW (
    -- * Memory allocation functions
    fftwMalloc,
    fftwFree,
    fftwFreePtr,
    fftwAllocReal,
    fftwAllocComplex,

    -- * FFT planning flags
    Direction(..),
    Flag(),
    -- ** Planning rigor flags
    fftwMeasure,
    fftwExhaustive,
    fftwPatient,
    fftwEstimate,
    fftwWisdomOnly,
    -- ** Algorithm restriction flags
    fftwDestroyInput,
    fftwUnaligned,
    fftwPreserveInput,

    -- * FFT planning
    FFTWPlan,
    planDFT1d,
    planDFTR2C1d,

    -- * FFT execution
    execute,
    executeDFT,
    executeDFTR2C
    ) where

import Foreign.C.Types
import Foreign.Ptr
import Data.Word
import Data.Complex
import Control.Monad
import Data.Bits
import Data.Monoid
import Data.Semigroup as Sem

#include <fftw3.h>

foreign import ccall unsafe "fftw_malloc"
    c_fftwMalloc :: CUInt -> IO (Ptr a)

-- | Like malloc, but ensures that the pointer obeys the alignment restrictions of FFTW (e.g. for SIMD acceleration). You probably want to use 'fftwAllocReal' or 'fftwAllocComplex' instead.
fftwMalloc :: Word32 -- ^ size
           -> IO (Ptr a)
fftwMalloc = c_fftwMalloc . fromIntegral

-- | Free a pointer returned by 'fftwMalloc', 'fftwAllocReal', or 'fftwAllocComplex'
foreign import ccall unsafe "fftw_free"
    fftwFree :: Ptr a -- ^ the pointer to be freed
             -> IO ()

-- | A function pointer to @fftwFree@.
foreign import ccall unsafe "&fftw_free"
    fftwFreePtr :: FunPtr (Ptr a -> IO ())

foreign import ccall unsafe "fftw_alloc_real"
    c_fftwAllocReal :: CUInt -> IO (Ptr CDouble)

-- | Allocates an array of Doubles. It ensures that the pointer obeys the alignment restrictions of FFTW (e.g. for SIMD acceleration).
fftwAllocReal :: Word32 -- ^ size
              -> IO (Ptr CDouble)
fftwAllocReal = c_fftwAllocReal . fromIntegral

foreign import ccall unsafe "fftw_alloc_complex"
    c_fftwAllocComplex :: CUInt -> IO (Ptr (Complex CDouble))

-- | Allocates an array of complex Doubles (i.e. the c type "double complex"). It ensures that the pointer obeys the alignment restrictions of FFTW (e.g. for SIMD acceleration).
fftwAllocComplex :: Word32 -- ^ size
                 -> IO (Ptr (Complex CDouble))
fftwAllocComplex = c_fftwAllocComplex . fromIntegral

-- | The direction of the transform: Forward for a normal transform, Backward for an inverse transform
data Direction = Forward
               | Backward

dirToInt :: Direction -> CInt
dirToInt Forward  = #const FFTW_FORWARD
dirToInt Backward = #const FFTW_BACKWARD

-- | FFTW planner flags. These flags affect the planning process. They can be combined using the 'Monoid' instance. See the FFTW flag documentation: <http://www.fftw.org/doc/Planner-Flags.html>.
newtype Flag = Flag {unFlag :: CUInt}

instance Sem.Semigroup Flag where
  (Flag x) <> (Flag y) = Flag (x .|. y)

instance Monoid Flag where
    mempty   = Flag 0
#if !(MIN_VERSION_base(4,11,0))
    mappend  = (Sem.<>)
#endif

fftwMeasure, fftwExhaustive, fftwPatient, fftwEstimate, fftwWisdomOnly :: Flag
fftwEstimate       = Flag #const FFTW_ESTIMATE
fftwMeasure        = Flag #const FFTW_MEASURE
fftwPatient        = Flag #const FFTW_PATIENT
fftwExhaustive     = Flag #const FFTW_EXHAUSTIVE
fftwWisdomOnly     = Flag #const FFTW_WISDOM_ONLY

fftwDestroyInput, fftwUnaligned, fftwPreserveInput :: Flag
fftwDestroyInput   = Flag #const FFTW_DESTROY_INPUT
fftwUnaligned      = Flag #const FFTW_UNALIGNED
fftwPreserveInput  = Flag #const FFTW_PRESERVE_INPUT

data CFFTWPlan

-- | A @FFTWPlan i o@ contains all of the information necessary to perform a transform from an input array of type @i@ to an output array of type @o@, including pointers to the input and output arrays.
newtype FFTWPlan i o = FFTWPlan (Ptr CFFTWPlan)

foreign import ccall unsafe "fftw_plan_dft_1d"
    c_planDFT1d :: CInt -> Ptr (Complex CDouble) -> Ptr (Complex CDouble) -> CInt -> CUInt -> IO (Ptr CFFTWPlan)

--This appears to be missing from the fft package on Hackage
-- | Create a plan for a 1 dimensional complex to complex DFT. The plan stores pointers to the input and output arrays, and these will be used if you 'execute' the plan in the future. They are required even if you intend to specify different input and output arrays in the future (i.e. using 'executeDFT')
planDFT1d :: Int                   -- ^ size
          -> Ptr (Complex CDouble) -- ^ input pointer
          -> Ptr (Complex CDouble) -- ^ output pointer
          -> Direction             -- ^ direction
          -> Flag                  -- ^ planner flags
          -> IO (FFTWPlan (Complex CDouble) (Complex CDouble))
planDFT1d n inp out sign flags = liftM FFTWPlan $ c_planDFT1d (fromIntegral n) inp out (dirToInt sign) (unFlag flags)

foreign import ccall unsafe "fftw_plan_dft_r2c_1d"
    c_planDFTR2C1d :: CInt -> Ptr CDouble -> Ptr (Complex CDouble) -> CUInt -> IO (Ptr CFFTWPlan)

--This appears to be missing from the fft package on Hackage
-- | Create a plan for a 1 dimensional real to complex DFT. The plan stores pointers to the input and output arrays, and these will be used if you 'execute' the plan in the future. They are required even if you intend to specify different input and output arrays in the future (i.e. using 'executeDFTR2C')
planDFTR2C1d :: Int                   -- ^ size
             -> Ptr CDouble           -- ^ input pointer
             -> Ptr (Complex CDouble) -- ^ output pointer
             -> Flag                  -- ^ planner flags
             -> IO (FFTWPlan CDouble (Complex CDouble))
planDFTR2C1d n inp out flags = liftM FFTWPlan $ c_planDFTR2C1d (fromIntegral n) inp out (unFlag flags)

foreign import ccall unsafe "fftw_execute"
    c_execute :: Ptr CFFTWPlan -> IO ()

-- | Execute a plan. Performs an FFT. The input and output arrays are stored within the plan so do not need to be given.
execute :: FFTWPlan i o -- ^ the plan to execute
        -> IO ()
execute (FFTWPlan p) = c_execute p

foreign import ccall unsafe "fftw_execute_dft"
    c_executeDFT :: Ptr CFFTWPlan -> Ptr (Complex CDouble) -> Ptr (Complex CDouble) -> IO ()

-- | Execute a complex to complex DFT but on different input and output arrays to those specified when the plan was created.
executeDFT :: FFTWPlan (Complex CDouble) (Complex CDouble) -- ^ the plan to execute
           -> Ptr (Complex CDouble)                        -- ^ input pointer
           -> Ptr (Complex CDouble)                        -- ^ output pointer
           -> IO ()
executeDFT (FFTWPlan p) inp out = c_executeDFT p inp out

foreign import ccall unsafe "fftw_execute_dft_r2c"
    c_executeDFTR2C :: Ptr CFFTWPlan -> Ptr CDouble -> Ptr (Complex CDouble) -> IO ()

-- | Execute a real to complex DFT but on different input and output arrays to those specified when the plan was created.
executeDFTR2C :: FFTWPlan CDouble (Complex CDouble) -- ^ the plan to execute
              -> Ptr CDouble                        -- ^ input pointer
              -> Ptr (Complex CDouble)              -- ^ output pointer
              -> IO ()
executeDFTR2C (FFTWPlan p) inp out = c_executeDFTR2C p inp out

ext: hsc

max_stars_repo_name: adamwalker/haskell-fftw-simple

max_issues_repo_name: adamwalker/haskell-fftw-simple

max_forks_repo_name: adamwalker/haskell-fftw-simple


text:
{-# LANGUAGE ScopedTypeVariables, TypeFamilies #-}

module Fuml.Supervised.KNN where

import Numeric.LinearAlgebra
import Data.List (sortBy)
import Data.Ord (comparing)
import Fuml.Core
import Fuml.Base.KNN
import Control.Monad.Identity

data KNN a = KNN {
   knnPoints :: [(Vector Double, a)]
  }

generalKNN :: Int -> (Vector Double -> Vector Double -> Double) -> ([(Double, a)] -> b) -> Supervisor Identity a (KNN a) b
generalKNN k dist vote  = Supervisor $ \_ theData ->
  return $ Predict (KNN theData) $ \v ->
   let ptDist (v1, y) = (dist v v1,y)
       nns = take k $ sortBy (comparing fst) $ map ptDist $ theData
   in vote nns

binaryKNN :: Int -> Supervisor Identity Bool (KNN Bool) Double
binaryKNN k = generalKNN k euclideanDistance weightedBoolVote

catKNN :: Eq a => Int -> Supervisor Identity a (KNN a) a
catKNN k = generalKNN k euclideanDistance weightedMajorityVote

regressKNN :: Int -> Supervisor Identity Double (KNN Double) Double
regressKNN k = generalKNN k euclideanDistance weightedAvgVote

max_stars_repo_name: ekalosak/open

max_issues_repo_name: ekalosak/open

max_forks_repo_name: ekalosak/open


text:
module HLearn.Metrics.Mahalanobis
    where

import Control.DeepSeq
import qualified Data.Vector.Generic as VG

import qualified Data.Strict as Strict
import Foreign.Storable
import Numeric.LinearAlgebra hiding ((<>),dim)
import qualified Numeric.LinearAlgebra as LA

import HLearn.Algebra

-------------------------------------------------------------------------------
-- data types

data Mahalanobis dp = Mahalanobis
    { rawdp :: !dp
    , moddp :: !dp
    }
    deriving (Read,Show,Eq,Ord)

instance NFData dp => NFData (Mahalanobis dp) where
    rnf m = deepseq (rawdp m)
          $ rnf (moddp m)

mkIdentity :: dp -> Mahalanobis dp
mkIdentity dp = Mahalanobis dp dp

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

dim :: VG.Vector dp r => Mahalanobis (dp r) -> Int
dim m = VG.length $ rawdp m

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

class
    ( Storable (Ring dp)
    , Element (Ring dp)
    , Container Vector (Ring dp)
    , LA.Product (Ring dp)
    , Field (Ring dp)
    , Ring dp ~ Double
    ) => MatrixField dp

instance
    ( Storable (Ring dp)
    , Element (Ring dp)
    , Container Vector (Ring dp)
    , LA.Product (Ring dp)
    , Field (Ring dp)
    , Ring dp ~ Double
    ) => MatrixField dp

class MkMahalanobis params where
    type MetricDatapoint params -- = Datapoint params
    mkMahalanobis :: params -> MetricDatapoint params -> Mahalanobis (MetricDatapoint params)

-------------------------------------------------------------------------------
-- algebra

instance HasRing dp => HasRing (Mahalanobis dp) where
    type Ring (Mahalanobis dp) = Ring dp

instance
    ( RealFrac r
    , Floating r
    , Ring (dp r) ~ r
    , HasRing (dp r)
    , VG.Vector dp r
    ) => MetricSpace (Mahalanobis (dp r)) where
    {-# INLINABLE distance #-}

    distance !m1 !m2 = {-# SCC distance #-} sqrt $ go 0 (dim m1-1)
        where
            go tot (-1) = tot
            go tot i = go (tot+(((rawdp m1) `VG.unsafeIndex` i)-((rawdp m2) `VG.unsafeIndex` i))
                              *(((moddp m1) `VG.unsafeIndex` i)-((moddp m2) `VG.unsafeIndex` i))) (i-1)

--     isFartherThanWithDistance !m1 !m2 !dist = {-# SCC isFartherThanWithDistance #-}
--         go 0 (dim m1-1)
--         where
--             dist2=dist*dist
--
--             go tot (-1) = Strict.Just $ sqrt tot
--             go tot i = if tot'>dist2
--                 then Strict.Nothing
--                 else go tot' (i-1)
--                 where
--                     tot' = tot+(((rawdp m1) `VG.unsafeIndex` i)-((rawdp m2) `VG.unsafeIndex` i))
--                               *(((rawdp m1) `VG.unsafeIndex` i)-((rawdp m2) `VG.unsafeIndex` i))

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
module Data.Gaussian where

import qualified Data.Complex as Complex
  -- Enables toComplex.
  -- Qualified to avoid clashing.

data Gaussian = (:+)
  {
    real :: !Integer,
    imaginary :: !Integer
  } deriving
    (Show, Read, Eq)
  -- The gaussian integer type.
  -- A gaussian integer is a complex number of the form a + bi where a and b are integers.

instance Num Gaussian
  where
    (a1 :+ b1) + (a2 :+ b2) = (a1 + a2) :+ (b1 + b2)
    (a1 :+ b1) - (a2 :+ b2) = (a1 - a2) :+ (b1 - b2)
    (a1 :+ b1) * (a2 :+ b2) = (a1 * a2 - b1 * b2) :+ (a1 * b2 + a2 * b1)
    abs (a :+ b) = ((a ^ two) + (b ^ two)) :+ 0
      where
        two = 2 :: Integer
          -- This prevents compiler warnings about defaulting types.
    signum = id
    negate (a :+ b) = negate a :+ negate b
    fromInteger x = x :+ 0
  -- The Num instance for Gaussian.
  -- We use the euclidean norm as abs: norm(a + bi) = a^2 + b^2

i :: Gaussian
i = 0 :+ 1
  -- The imaginary unit, equal to sqrt(-1).

magnitude :: Floating c => Gaussian -> c
magnitude (a :+ b) = sqrt $ (fromInteger a ^ two) + (fromInteger b ^ two)
  where
    two = 2 :: Integer
      -- This prevents compiler warnings about defaulting types.
  -- Equivalent to sqrt . fromInteger . real . abs

norm :: Gaussian -> Integer
norm (a :+ b) = (a ^ two) + (b ^ two)
  where
    two = 2 :: Integer
      -- This prevents compiler warnings about defaulting types.
  -- Equivalent to real . abs

conjugate :: Gaussian -> Gaussian
conjugate (a :+ b) = (a :+ negate b)
  -- The complex conjugate.
  -- Effectively flips about the real axis.

rotate_cw :: Gaussian -> Gaussian
rotate_cw (a :+ b) = (b :+ negate a)
  -- Rotate a gaussian integer clockwise.
  -- Equivalent to multiplying by i.

rotate_ccw :: Gaussian -> Gaussian
rotate_ccw (a :+ b) = (negate b :+ a)
  -- Rotate a gaussian integer counterclockwise.
  -- Equivalent to multiplying by -i.

flip_x :: Gaussian -> Gaussian
flip_x (a :+ b) = (negate a :+ b)
  -- Flips on imaginary axis.
  -- Analogous to conjugate.

flip_y :: Gaussian -> Gaussian
flip_y = conjugate
  -- Alias for conjugate.
  -- See conjugate.

swap_x_y :: Gaussian -> Gaussian
swap_x_y (a :+ b) = (b :+ a)
  -- Utility function, swaps real and imaginary axes.

flip_diag_2 :: Gaussian -> Gaussian
flip_diag_2 (a :+ b) = (negate b :+ negate a)
  -- Utility function, reflects on x = -y
  -- Equivalent to negate . swap_x_y

toComplex :: Num p => Gaussian -> Complex.Complex p
toComplex (a :+ b) = (fromInteger a Complex.:+ fromInteger b)
  -- Converts a Gaussian to any complex number: Num p => Complex p
  -- Fully qualified operators to avoid clashing.

max_stars_repo_name: schuelermine/gaussian

max_issues_repo_name: schuelermine/gaussian

max_forks_repo_name: schuelermine/gaussian


text:
-- This file is part of Quipper. Copyright (C) 2011-2014. Please see the
-- file COPYRIGHT for a list of authors, copyright holders, licensing,
-- and other details. All rights reserved.
--
-- ======================================================================

{-# LANGUAGE NoMonomorphismRestriction #-}

-- | This module contains an implementation of the oracle and its
-- automatic lifting to quantum circuits using Template Haskell.
module Algorithms.QLS.TemplateOracle where

import Data.Complex
import Libraries.Auxiliary

import Control.Monad

import QuipperLib.Arith hiding (template_symb_plus_)
import Algorithms.QLS.Utils
import Algorithms.QLS.QDouble
import Algorithms.QLS.RealFunc
import Algorithms.QLS.QSignedInt
import Algorithms.QLS.CircLiftingImport

import Quipper
import Quipper.CircLifting

-- | Lifted version of @'any'@ (local version).
build_circuit
local_any :: (a -> Bool) -> [a] -> Bool
local_any f l = case l of
            []    -> False
            (h:t) -> if (f h) then True else local_any f t

-- | Auxiliary function.
build_circuit
itoxy :: Int -> Int -> Int -> (Int, Int)
itoxy i nx ny =
   if (i <= 0) then (0,0)
   else if (i > (nx-1)*ny + nx*(ny-1)) then (0,0)
   else ((mod (i - 1) (2*nx - 1)) + 1, ceiling ((fromIntegral i)/(2.0*(fromIntegral nx)-1.0)))

-- | The function sin /x/ \/ /x/.
build_circuit
sinc :: Double -> Double
sinc x = if (x /= 0.0) then (sin x) / x else 1.0

-- | Auxiliary function.
build_circuit
edgetoxy :: Int -> Int -> Int -> (Double, Double)
edgetoxy e nx ny =
  let (ex,ey) = itoxy e nx ny in
     if (ex < nx) then ((fromIntegral ex) + 0.5, fromIntegral ey)
     else (fromIntegral (ex - nx + 1), (fromIntegral ey) + 0.5)

-- | Auxiliary function. The inputs are:
--
-- * /y1/ :: 'Int' - Global edge index of row index of desired matrix element;
--
-- * /y2/ :: 'Int' - Global edge index of column index of desired matrix element;
--
-- * /nx/ :: 'Int' - Number of vertices left to right;
--
-- * /ny/ :: 'Int' - Number of vertices top to bottom;
--
-- * /lx/ :: 'Double' - Length of horizontal edges (distance between vertices in /x/ direction);
--
-- * /ly/ :: 'Double' - Length of vertical edges (distance between vertices in /y/ direction);
--
-- * /k/ :: 'Double' - Plane wave wavenumber.
--
-- The output is the matrix element /A/(/y1/, /y2/).
build_circuit
calcmatrixelement :: --
	-- Inputs
	Int 	-- int y1 - Global edge index of row index of desired matrix element
	-> Int 	-- int y2 - Global edge index of column index of desired matrix element
	-> Int 	-- int nx - Number of vertices left to right
	-> Int	-- int ny - Number of vertices top to bottom
	-> Double -- float lx - Length of horizontal edges (distance between vertices in x direction)
	-> Double -- float ly - Length of vertical edges (distance between vertices in y direction)
	-> Double -- float k - Plane wave wavenumber
	-- Outputs
	-> Complex Double -- float A - Matrix element A(y1,y2)
calcmatrixelement y1 y2 nx ny lx ly k =
	let (xg1, yg1) = itoxy y1 nx ny in
        let (xg2, yg2) = itoxy y2 nx ny in
        let b = if ( (y1==y2) && (xg1 >= nx) ) then ly/lx - k*k*lx*ly/3.0
                -- B11 and B22
                else if ( (y1==y2) && (xg1<nx) ) then lx/ly - k*k*lx*ly/3.0
                -- B12 and B21
                else if ( (abs(yg1-yg2) == 1) && (abs(xg1-xg2) == 0) && (xg1<nx) ) then -lx/ly + k*k*lx*ly/12.0
                -- B34 and B43
                else if ( (abs(yg1-yg2)==0) && (abs(xg1-xg2) == 1) && (xg1>=nx) ) then -ly/lx + k*k*lx*ly/12.0
                -- B13
                else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx+1)) && (xg2>=nx) ) then -1.0
                -- B31
                else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx+1)) && (xg1>=nx) ) then -1.0
                -- B14
                else if ( (yg1==(yg2+1)) && (xg1==(xg2-nx)) && (xg1<nx) ) then 1.0
                -- B41
                else if ( (yg2==(yg1+1)) && (xg2==(xg1-nx)) && (xg2<nx) ) then 1.0
                -- B42
                else if ( (yg1==yg2) && (xg1==(xg2+nx)) && (xg1>=nx) ) then -1.0
                -- B24
                else if ( (yg2==yg1) && (xg2==(xg1+nx)) && (xg2>=nx) ) then -1.0
                -- B32
                else if ( (yg1==yg2) && (xg1==(xg2+nx-1)) && (xg2<nx) ) then 1.0
                -- B23
                else if ( (yg2==yg1) && (xg2==(xg1+nx-1)) && (xg1<nx) ) then 1.0
                else -1.0
        in
        let c = if ( (y1==y2) && ( (xg1==nx) || (xg1==(2*nx-1)) ) ) then 0.0 :+ (k*ly)
                else if  ( (y1==y2) && ( ((yg1==1) && (xg1<nx)) || ((yg1==ny) && (xg1<nx)) ) ) then 0.0 :+ (k*lx)
                else 0.0 :+ 0.0
	in (b :+ 0.0) + c

-- | Auxiliary function.
build_circuit
get_edges l = case l of
                [] -> []
                (x:tt) -> case tt of
                   [] -> []
                   (y:t) -> (x,y):(get_edges (y:t))

-- | Auxiliary function.
build_circuit
checkedge :: Int -> [(Double,Double)] -> Int -> Int -> Bool
checkedge e scatteringnodes nx ny =
     let (xi,yi) = edgetoxy e nx ny in
     let test_elt ((x1,y1),(x2, y2)) = xi >= x1 && yi >= y1 && xi <= x2 && yi <= y2 in
     let half = take_half scatteringnodes in
     let test_list = local_any test_elt (get_edges half) in (not test_list)

-- | Oracle /r/.
build_circuit
calcRweights :: Int -> Int -> Int -> Double -> Double -> Double -> Double -> Double -> Complex Double
calcRweights y nx ny lx ly k theta phi =
     let (xc',yc') = edgetoxy y nx ny in
     let xc = (xc'-1.0)*lx - ((fromIntegral nx)-1.0)*lx/2.0 in
     let yc = (yc'-1.0)*ly - ((fromIntegral ny)-1.0)*ly/2.0 in
     let (xg,yg) = itoxy y nx ny in

     if (xg == nx) then

         let i = (mkPolar ly (k*xc*(cos phi)))*
                 (mkPolar 1.0 (k*yc*(sin phi)))*
                 ((sinc (k*ly*(sin phi)/2.0)) :+ 0.0) in

         let r = ( cos(phi) :+ k*lx )*((cos (theta - phi))/lx :+ 0.0) in i * r

     else if (xg==2*nx-1) then

         let i = (mkPolar ly (k*xc*cos(phi)))*
                 (mkPolar 1.0 (k*yc*sin(phi)))*
                 ((sinc (k*ly*sin(phi)/2.0)) :+ 0.0) in

         let r = ( cos(phi) :+ (- k*lx))*((cos (theta - phi))/lx :+ 0.0) in i * r

     else if ( (yg==1) && (xg<nx) ) then

         let i = (mkPolar lx (k*yc*sin(phi)))*
                 (mkPolar 1.0 (k*xc*cos(phi)))*
                 ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

         let r = ( (- sin phi) :+ k*ly )*((cos(theta - phi))/ly :+ 0.0) in i * r

     else if ( (yg==ny) && (xg<nx) ) then

         let i = (mkPolar lx (k*yc*sin(phi)))*
                 (mkPolar 1.0 (k*xc*cos(phi)))*
                 ((sinc (k*lx*(cos phi)/2.0)) :+ 0.0) in

         let r = ( (- sin phi) :+ (- k*ly) )*((cos(theta - phi)/ly) :+ 0.0) in i * r

     else 0.0 :+ 0.0

-- | Auxiliary function for oracle /A/.
build_circuit
convertband :: Int -> Int -> Int -> Int -> Int
convertband y b nx ny =
 let nedges = (nx - 1)*ny + nx*(ny - 1) in
 let (ex,ey) = itoxy y nx ny in
 let x = if ( (ex < nx) && (ey /= 1) ) then
           case b of
                1 -> y-2*nx+1
                2 -> y-nx
                3 -> y-nx+1
                5 -> y
                7 -> y+nx-1
                8 -> y+nx
                9 -> y+2*nx-1
                _ -> -1
         else if ( (ex < nx) && (ey == 1) ) then
           case b of
                5 -> y
                7 -> y+nx-1
                8 -> y+nx
                9 -> y+2*nx-1
                _ -> -1
         else if ( (ex >= nx) && (ex /= nx) && (ex /= 2*nx-1) ) then
           case b of
                    2 -> y-nx
                    3 -> y-nx+1
                    4 -> y-1
                    5 -> y
                    6 -> y+1
                    7 -> y+nx-1
                    8 -> y+nx
                    _ -> -1
         else if ( (ex >= nx) && (ex == nx) ) then
                 case b of
                    3 -> y-nx+1
                    5 -> y
                    6 -> y+1
                    8 -> y+nx
                    _ -> -1
         else if ( (ex >= nx) && (ex == 2*nx-1) ) then
                case b of
                    2 -> y-nx
                    4 -> y-1
                    5 -> y
                    7 -> y+nx-1
                    _ -> -1
         else -1
 in if ( (x < 1) || (x > nedges) ) then -1 else x

-- | Oracle /A/. It is equivalent to the Matlab function
-- 'getBandNodeValues'.
--
-- 'getNodeValuesMoreOutputs v b ...'  outputs the node of the edge
-- connected to vertex v in band b, and a real number parameterized by
-- the 'BoolParam' parameter: the magnitude (PFalse) or the phase
-- (PTrue) of the complex value at the corresponding place in the matrix A.
build_circuit
getNodeValuesMoreOutputs ::
    Int -> Int -> Int -> Int -> [(Double,Double)] -> Double -> Double -> Double -> BoolParam
    -> Int -> (Int, Double)
getNodeValuesMoreOutputs v' b' nx ny scatteringnodes lx ly k argflag maxConnectivity =
   let maxC = getIntFromParam maxConnectivity in
   let b = getIntFromParam b' in
   let nedges = (nx - 1)*ny + nx*(ny - 1) in
   let flag = v' <= nedges in
   let v = (if flag then v' else (v' - nedges)) in
   let nodeDefault = (if flag then v + nedges else v) in
   let valueDefault = 0.0 in
   let indicesDefault = (-1, -1) in
   let isvalid = checkedge v scatteringnodes nx ny
   in
   if ( (not isvalid) && b == 5 ) then

     let indices = (if flag then (v, v + nedges) else (v + nedges, v)) in
     case argflag of
        PTrue  -> (nodeDefault, valueDefault)
        PFalse -> (nodeDefault, 1.0)

   else if ( (not isvalid) && b /= 5) then (nodeDefault, valueDefault)

   else if ((b > maxC + 2) || (b <= 0)) then (nodeDefault, valueDefault)

   else

   let x = convertband v b' nx ny in

   if (x == -1) then (nodeDefault, valueDefault)
   else

   let isvalid = checkedge x scatteringnodes nx ny in

   if isvalid then

     let ax = calcmatrixelement v x nx ny lx ly k in
     let (node, indices) = if flag then (x + nedges, (v, x + nedges))
                           else (x, (v+nedges,x)) in
     let value = case argflag of
                   PTrue -> atan2  (imagPart ax) (realPart ax)
                   PFalse -> magnitude ax
     in (node, value)

   else
     (nodeDefault, valueDefault)

-- | Auxiliary function for oracle /b/. The inputs are:
--
-- * /y/ :: 'Int' - Global edge index.  Note this is the unmarked /y/
-- coordinate, i.e. the coordinate without scattering regions removed;
--
-- * /nx/ :: 'Int' - Number of vertices left to right;
--
-- * /ny/ :: 'Int' - Number of vertices top to bottom;
--
-- * /lx/ :: 'Double' - Length of horizontal edges (distance between vertices in /x/ direction);
--
-- * /ly/ :: 'Double' - Length of vertical edges (distance between vertices in /y/ direction);
--
-- * /k/ :: 'Double' - Plane wave wavenumber;
--
-- * θ :: 'Double' - Direction of wave propagation;
--
-- * /E0/ :: 'Double' - Magnitude of incident plane wave.
--
-- The output is the magnitude of the electric field on edge /y/.
build_circuit
calcincidentfield :: --
	-- Inputs
		Int -- int y - Global edge index.  Note this is the unmarked y coordinate,
			-- i.e. the coordinate without scattering regions removed.
	-> Int -- int nx - Number of vertices left to right
	-> Int --  int ny - Number of vertices top to bottom
	-> Double -- float lx - Length of horizontal edges (distance between vertices in x direction)
	-> Double -- float ly - Length of vertical edges (distance between vertices in y direction)
	-> Double -- float k - Plane wave wavenumber
	-> Double -- float theta - Direction of wave propagation
	-> Double -- float E0 - Magnitude of incident plane wave
	-- Outputs
	-> Complex Double -- complex float e - Magnitude of electric field on edge y
calcincidentfield y nx ny lx ly k theta e0 =
        let (xg, yg) = itoxy y nx ny in
        --Determine whether edge is horizontal or vertical
        let isvertical = xg >= nx in
        let (xvalueTmp, yvalueTmp) = edgetoxy y nx ny in
        let xvalue = xvalueTmp * lx in
        let yvalue = yvalueTmp * ly in
        -- Convert x and y edge coordinates to x and y values and caluculate field
        if isvertical then
          mkPolar (-cos(theta)*e0) ( -k*(xvalue*cos(theta)+yvalue*sin(theta)))
        else
          mkPolar (sin(theta)*e0) ( -k*(xvalue*cos(theta)+yvalue*sin(theta)))

-- | Auxiliary function for oracle /b/.
build_circuit
getconnection :: Int -> Int -> Int -> Int -> Int -> Int
getconnection y i' nx ny maxConnectivity =
  let i = getIntFromParam i' in
  let maxC = getIntFromParam maxConnectivity in
  let (ex,ey) = itoxy y nx ny in
  let x = if ( (ex < nx) && (ey /= 1) ) then
            case i' of
                1 -> y-2*nx+1
                2 -> y-nx
                3 -> y-nx+1
                4 -> y
                5 -> y+nx-1
                6 -> y+nx
                7 -> y+2*nx-1
                _ -> -1
          else if ( (ex < nx) && (ey == 1) ) then
            case i' of
                1 -> y
                2 -> y+nx-1
                3 -> y+nx
                4 -> y+2*nx-1
                _ -> -1
          else if ( (ex >= nx) && (ex /= nx) && (ex /= 2*nx-1) ) then
             case i' of
                1 -> y-nx
                2 -> y-nx+1
                3 -> y-1
                4 -> y
                5 -> y+1
                6 -> y+nx-1
                7 -> y+nx
                _ -> -1
          else if ( (ex >= nx) && (ex == nx) ) then
             case i' of
                1 -> y-nx+1
                2 -> y
                3 -> y+1
                4 -> y+nx
                _ -> -1
          else if ( (ex >= nx) && (ex == 2*nx-1) ) then
             case i' of
                1 -> y-nx
                2 -> y-1
                3 -> y
                4 -> y+nx-1
                _ -> -1
          else -1
  in
  if (i > maxC) then -1
  else if (x > nx*(ny-1)+ny*(nx-1)) then -1
  else x

-- | Auxiliary function to @'template_paramZero'@.
build_circuit
local_loop_with_index_aux :: Int -> Int -> t -> (Int -> t -> t) -> t
local_loop_with_index_aux i n x f =
   case paramMinus n i of
     0 -> x
     _ -> local_loop_with_index_aux (paramSucc i) n (f i x) f

-- | Local version of @'loop_with_index'@, for lifting.
build_circuit
local_loop_with_index :: Int -> t -> (Int -> t -> t) -> t
local_loop_with_index n x f = local_loop_with_index_aux paramZero n x f

-- | Oracle /b/.
build_circuit
getKnownWeights :: Int -> Int -> Int -> [(Double,Double)] -> Double -> Double -> Double -> Double -> Double -> Int -> Complex Double
getKnownWeights y nx ny scatteringnodes lx ly k theta e0 maxConnectivity =
   let makeConnections i connections = let x = getconnection y (paramSucc i) nx ny maxConnectivity in
                                       let t = not $ checkedge x scatteringnodes nx ny in
                                       (x,t):connections
   in
   let calcTang b (c,t) = if t
                          then let matElt = calcmatrixelement y c nx ny lx ly k in
                               let incField = calcincidentfield c nx ny lx ly k theta e0
                               in b - matElt * incField
                          else b
   in
   let connections = local_loop_with_index maxConnectivity [] makeConnections in
   if (not $ checkedge y scatteringnodes nx ny) then 0.0 :+ 0.0
   else foldl calcTang (0.0 :+ 0.0) connections

----------------------------------------------------------------------
----------------------------------------------------------------------
-- Testing functions

test_template_sinc = do
      f <- template_sinc
      r <- qinit (0 :: FDouble)
      f r
      return ()

test_template_itoxy = do
      f <- template_itoxy
      x <- qinit (0 :: FSignedInt)
      y <- qinit (0 :: FSignedInt)
      z <- qinit (0 :: FSignedInt)
      g <- f x
      h <- g y
      k <- h z
      return ()

test_template_edgetoxy = do
      f <- template_edgetoxy
      x <- qinit (0 :: FSignedInt)
      y <- qinit (0 :: FSignedInt)
      z <- qinit (0 :: FSignedInt)
      g <- f x
      h <- g y
      k <- h z
      return ()

test_template_calcRweights = do
      f <- template_calcRweights
      n1 <- qinit (0 :: FSignedInt)
      n2 <- qinit (0 :: FSignedInt)
      n3 <- qinit (0 :: FSignedInt)
      x1 <- qinit (0 :: FDouble)
      x2 <- qinit (0 :: FDouble)
      x3 <- qinit (0 :: FDouble)
      x4 <- qinit (0 :: FDouble)
      x5 <- qinit (0 :: FDouble)
      f1 <- f n1
      f2 <- f1 n2
      f3 <- f2 n3
      g1 <- f3 x1
      g2 <- g1 x2
      g3 <- g2 x3
      g4 <- g3 x4
      g5 <- g4 x5
      return ()

test_template_calcincidentfield = do
   y' <- qinit (0 :: FSignedInt)
   nx' <- qinit (0 :: FSignedInt)
   ny' <- qinit (0 :: FSignedInt)
   lx' <- qinit (0 :: FDouble)
   ly' <- qinit (0 :: FDouble)
   k' <- qinit (0 :: FDouble)
   theta' <- qinit (0 :: FDouble)
   e0' <- qinit (0 :: FDouble)
   f <- template_calcincidentfield
   f1 <- f y'
   f2 <- f1 nx'
   f3 <- f2 ny'
   f4 <- f3 lx'
   f5 <- f4 ly'
   f6 <- f5 k'
   f7 <- f6 theta'
   f8 <- f7 e0'
   return ()

test_template_calcmatrixelement = do
   y1 <- qinit (0 :: FSignedInt)
   y2 <- qinit (0 :: FSignedInt)
   nx <- qinit (0 :: FSignedInt)
   ny <- qinit (0 :: FSignedInt)
   lx <- qinit (0 :: FDouble)
   ly <- qinit (0 :: FDouble)
   k  <- qinit (0 :: FDouble)
   f <- template_calcmatrixelement
   f1 <- f y1
   f2 <- f1 y2
   f3 <- f2 nx
   f4 <- f3 ny
   f5 <- f4 lx
   f6 <- f5 ly
   f7 <- f6 k
   return ()

test_template_getconnection = do
   y <- qinit (0 :: FSignedInt)
   nx <- qinit (0 :: FSignedInt)
   ny <- qinit (0 :: FSignedInt)
   f <- template_getconnection
   f1 <- f y
   f2 <- f1 6
   f3 <- f2 nx
   f4 <- f3 ny
   f5 <- f4 7
   return ()

test_template_checkedge = do
   y <- qinit (0 :: FSignedInt)
   nx <- qinit (0 :: FSignedInt)
   ny <- qinit (0 :: FSignedInt)
   s <- qinit [(0 :: FDouble,0 :: FDouble),(0 :: FDouble,0 :: FDouble)]
   f <- template_checkedge
   f1 <- f y
   f2 <- f1 s
   f3 <- f2 nx
   f4 <- f3 ny
   return ()

max_stars_repo_name: fritzo/quipper

max_issues_repo_name: fritzo/quipper

max_forks_repo_name: fritzo/quipper


text:
{-|
Copyright  :  (C) 2013-2016, University of Twente,
                  2016-2017, Myrtle Software Ltd
License    :  BSD2 (see the file LICENSE)
Maintainer :  Christiaan Baaij <christiaan.baaij@gmail.com>
-}

{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}

{-# LANGUAGE Trustworthy #-}

{-# OPTIONS_GHC -fplugin=GHC.TypeLits.Extra.Solver #-}
{-# OPTIONS_GHC -fplugin=GHC.TypeLits.Normalise #-}
{-# OPTIONS_GHC -fplugin=GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_HADDOCK show-extensions #-}

#include "MachDeps.h"

module Clash.Class.BitPack
  ( BitPack (..)
  , bitCoerce
  , bitCoerceMap
  , boolToBV
  , boolToBit
  , bitToBool
  , packXWith

  -- * Internals
  , GBitPack(..)
  )
where

import Control.Exception              (catch, evaluate)
import Data.Binary.IEEE754            (doubleToWord, floatToWord, wordToDouble,
                                       wordToFloat)

#if MIN_VERSION_base(4,12,0)
import Data.Complex                   (Complex)
import Data.Ord                       (Down)
#endif

import Data.Int
import Data.Word
import Foreign.C.Types                (CUShort)
import GHC.TypeLits                   (KnownNat, Nat, type (+), type (-))
import Numeric.Half                   (Half (..))
import GHC.Generics
import GHC.TypeLits.Extra             (CLog, Max)
import Prelude                        hiding (map)
import System.IO.Unsafe               (unsafeDupablePerformIO)

import Clash.Promoted.Nat             (SNat(..), snatToNum)
import Clash.Class.BitPack.Internal   (deriveBitPackTuples)
import Clash.Class.Resize             (zeroExtend, resize)
import Clash.Sized.BitVector          (Bit, BitVector, (++#))
import Clash.Sized.Internal.BitVector
  (pack#, split#, checkUnpackUndef, undefined#, unpack#, unsafeToNatural)
import Clash.XException

{- $setup
>>> :set -XDataKinds
>>> import Clash.Prelude
-}

-- | Convert to and from a 'BitVector'
class KnownNat (BitSize a) => BitPack a where
  -- | Number of 'Clash.Sized.BitVector.Bit's needed to represents elements
  -- of type @a@
  --
  -- Can be derived using `GHC.Generics`:
  --
  -- > import Clash.Prelude
  -- > import GHC.Generics
  -- >
  -- > data MyProductType = MyProductType { a :: Int, b :: Bool }
  -- >   deriving (Generic, BitPack)
  type BitSize a :: Nat
  type BitSize a = (CLog 2 (GConstructorCount (Rep a))) + (GFieldSize (Rep a))
  -- | Convert element of type @a@ to a 'BitVector'
  --
  -- >>> pack (-5 :: Signed 6)
  -- 11_1011
  pack   :: a -> BitVector (BitSize a)
  default pack
    :: ( Generic a
       , GBitPack (Rep a)
       , KnownNat (BitSize a)
       , KnownNat constrSize
       , KnownNat fieldSize
       , constrSize ~ CLog 2 (GConstructorCount (Rep a))
       , fieldSize ~ GFieldSize (Rep a)
       , (constrSize + fieldSize) ~ BitSize a
       )
    => a -> BitVector (BitSize a)
  pack = packXWith go
   where
    go a = resize (pack sc) ++# packedFields
     where
      (sc, packedFields) = gPackFields 0 (from a)

  -- | Convert a 'BitVector' to an element of type @a@
  --
  -- >>> pack (-5 :: Signed 6)
  -- 11_1011
  -- >>> let x = pack (-5 :: Signed 6)
  -- >>> unpack x :: Unsigned 6
  -- 59
  -- >>> pack (59 :: Unsigned 6)
  -- 11_1011
  unpack :: BitVector (BitSize a) -> a
  default unpack
    :: ( Generic a
       , GBitPack (Rep a)
       , KnownNat constrSize
       , KnownNat fieldSize
       , constrSize ~ CLog 2 (GConstructorCount (Rep a))
       , fieldSize ~ GFieldSize (Rep a)
       , (constrSize + fieldSize) ~ BitSize a
       )
    => BitVector (BitSize a) -> a
  unpack b =
    to (gUnpack sc 0 bFields)
   where
    (checkUnpackUndef unpack . resize -> sc, bFields) = split# b

packXWith
  :: KnownNat n
  => (a -> BitVector n)
  -> a
  -> BitVector n
packXWith f x =
  unsafeDupablePerformIO (catch (f <$> evaluate x)
                                (\(XException _) -> return undefined#))
{-# NOINLINE packXWith #-}

{-# INLINE[1] bitCoerce #-}
-- | Coerce a value from one type to another through its bit representation.
--
-- >>> pack (-5 :: Signed 6)
-- 11_1011
-- >>> bitCoerce (-5 :: Signed 6) :: Unsigned 6
-- 59
-- >>> pack (59 :: Unsigned 6)
-- 11_1011
bitCoerce
  :: (BitPack a, BitPack b, BitSize a ~ BitSize b)
  => a
  -> b
bitCoerce = unpack . pack

-- | Map a value by first coercing to another type through its bit representation.
--
-- >>> pack (-5 :: Signed 32)
-- 1111_1111_1111_1111_1111_1111_1111_1011
-- >>> bitCoerceMap @(Vec 4 (BitVector 8)) (replace 1 0) (-5 :: Signed 32)
-- -16711685
-- >>> pack (-16711685 :: Signed 32)
-- 1111_1111_0000_0000_1111_1111_1111_1011
bitCoerceMap
  :: forall a b . (BitPack a, BitPack b, BitSize a ~ BitSize b)
  => (a -> a)
  -> b
  -> b
bitCoerceMap f = bitCoerce . f . bitCoerce

instance BitPack Bool where
  type BitSize Bool = 1
  pack   = let go b = if b then 1 else 0 in packXWith go
  unpack = checkUnpackUndef $ \bv -> if bv == 1 then True else False

instance KnownNat n => BitPack (BitVector n) where
  type BitSize (BitVector n) = n
  pack     = packXWith id
  unpack v = v

instance BitPack Bit where
  type BitSize Bit = 1
  pack   = packXWith pack#
  unpack = unpack#

instance BitPack Int where
  type BitSize Int = WORD_SIZE_IN_BITS
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Int8 where
  type BitSize Int8 = 8
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Int16 where
  type BitSize Int16 = 16
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Int32 where
  type BitSize Int32 = 32
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Int64 where
  type BitSize Int64 = 64
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Word where
  type BitSize Word = WORD_SIZE_IN_BITS
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Word8 where
  type BitSize Word8 = 8
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Word16 where
  type BitSize Word16 = 16
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Word32 where
  type BitSize Word32 = 32
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Word64 where
  type BitSize Word64 = 64
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Float where
  type BitSize Float = 32
  pack   = packXWith packFloat#
  unpack = checkUnpackUndef unpackFloat#

packFloat# :: Float -> BitVector 32
packFloat# = fromIntegral . floatToWord
{-# NOINLINE packFloat# #-}

unpackFloat# :: BitVector 32 -> Float
unpackFloat# (unsafeToNatural -> w) = wordToFloat (fromIntegral w)
{-# NOINLINE unpackFloat# #-}

instance BitPack Double where
  type BitSize Double = 64
  pack   = packXWith packDouble#
  unpack = checkUnpackUndef unpackDouble#

packDouble# :: Double -> BitVector 64
packDouble# = fromIntegral . doubleToWord
{-# NOINLINE packDouble# #-}

unpackDouble# :: BitVector 64 -> Double
unpackDouble# (unsafeToNatural -> w) = wordToDouble (fromIntegral w)
{-# NOINLINE unpackDouble# #-}

instance BitPack CUShort where
  type BitSize CUShort = 16
  pack   = packXWith fromIntegral
  unpack = checkUnpackUndef fromIntegral

instance BitPack Half where
  type BitSize Half = 16
  pack (Half x) = pack x
  unpack        = checkUnpackUndef $ \x -> Half (unpack x)

instance BitPack () where
  type BitSize () = 0
  pack   _ = minBound
  unpack _ = ()

instance (BitPack a, BitPack b) =>
    BitPack (a,b) where
  type BitSize (a,b) = BitSize a + BitSize b
  pack = let go (a,b) = pack a ++# pack b in packXWith go
  unpack ab  = let (a,b) = split# ab in (unpack a, unpack b)

class GBitPack f where
  -- | Size of fields. If multiple constructors exist, this is the maximum of
  -- the sum of each of the constructors fields.
  type GFieldSize f :: Nat

  -- | Number of constructors this type has. Indirectly indicates how many bits
  -- are needed to represent the constructor.
  type GConstructorCount f :: Nat

  -- | Pack fields of a type. Caller should pack and prepend the constructor bits.
  gPackFields
    :: Int
    -- ^ Current constructor
    -> f a
    -- ^ Data to pack
    -> (Int, BitVector (GFieldSize f))
    -- ^ (Constructor number, Packed fields)

  -- | Unpack whole type.
  gUnpack
    :: Int
    -- ^ Construct with constructor /n/
    -> Int
    -- ^ Current constructor
    -> BitVector (GFieldSize f)
    -- ^ BitVector containing fields
    -> f a
    -- ^ Unpacked result

instance GBitPack a => GBitPack (M1 m d a) where
  type GFieldSize (M1 m d a) = GFieldSize a
  type GConstructorCount (M1 m d a) = GConstructorCount a

  gPackFields cc (M1 m1) = gPackFields cc m1
  gUnpack c cc b = M1 (gUnpack c cc b)

instance ( KnownNat (GFieldSize g)
         , KnownNat (GFieldSize f)
         , KnownNat (GConstructorCount f)
         , GBitPack f
         , GBitPack g
         ) => GBitPack (f :+: g) where
  type GFieldSize (f :+: g) = Max (GFieldSize f) (GFieldSize g)
  type GConstructorCount (f :+: g) = GConstructorCount f + GConstructorCount g

  gPackFields cc (L1 l) =
    let (sc, packed) = gPackFields cc l in
    let padding = undefined# :: BitVector (Max (GFieldSize f) (GFieldSize g) - GFieldSize f) in
    (sc, packed ++# padding)
  gPackFields cc (R1 r) =
    let cLeft = snatToNum (SNat @(GConstructorCount f)) in
    let (sc, packed) = gPackFields (cc + cLeft) r in
    let padding = undefined# :: BitVector (Max (GFieldSize f) (GFieldSize g) - GFieldSize g) in
    (sc, packed ++# padding)

  gUnpack c cc b =
    let cLeft = snatToNum (SNat @(GConstructorCount f)) in
    if c < cc + cLeft then
      L1 (gUnpack c cc f)
    else
      R1 (gUnpack c (cc + cLeft) g)

   where
    -- It's a thing of beauty, if I may say so myself!
    (f, _ :: BitVector (Max (GFieldSize f) (GFieldSize g) - GFieldSize f)) = split# b
    (g, _ :: BitVector (Max (GFieldSize f) (GFieldSize g) - GFieldSize g)) = split# b

instance (KnownNat (GFieldSize g), KnownNat (GFieldSize f), GBitPack f, GBitPack g) => GBitPack (f :*: g) where
  type GFieldSize (f :*: g) = GFieldSize f + GFieldSize g
  type GConstructorCount (f :*: g) = 1

  gPackFields cc fg =
    (cc, packXWith go fg)
   where
    go (l0 :*: r0) =
      let (_, l1) = gPackFields cc l0 in
      let (_, r1) = gPackFields cc r0 in
      l1 ++# r1

  gUnpack c cc b =
    gUnpack c cc front :*: gUnpack c cc back
   where
    (front, back) = split# b

instance BitPack c => GBitPack (K1 i c) where
  type GFieldSize (K1 i c) = BitSize c
  type GConstructorCount (K1 i c)  = 1

  gPackFields cc (K1 i) = (cc, pack i)
  gUnpack _c _cc b      = K1 (unpack b)

instance GBitPack U1 where
  type GFieldSize U1 = 0
  type GConstructorCount U1 = 1

  gPackFields cc U1 = (cc, 0)
  gUnpack _c _cc _b = U1

-- Instances derived using Generic
instance ( BitPack a
         , BitPack b
         ) => BitPack (Either a b)

instance BitPack a => BitPack (Maybe a)

#if MIN_VERSION_base(4,12,0)
instance BitPack a => BitPack (Complex a)
instance BitPack a => BitPack (Down a)
#endif

-- | Zero-extend a 'Bool'ean value to a 'BitVector' of the appropriate size.
--
-- >>> boolToBV True :: BitVector 6
-- 00_0001
-- >>> boolToBV False :: BitVector 6
-- 00_0000
boolToBV :: KnownNat n => Bool -> BitVector (n + 1)
boolToBV = zeroExtend . pack

-- | Convert a Bool to a Bit
boolToBit :: Bool -> Bit
boolToBit = bitCoerce

-- | Convert a Bool to a Bit
bitToBool :: Bit -> Bool
bitToBool = bitCoerce

-- Derive the BitPack instance for tuples of size 3 to 62
deriveBitPackTuples ''BitPack ''BitSize 'pack 'unpack

max_stars_repo_name: gergoerdi/clash-compiler

max_issues_repo_name: gergoerdi/clash-compiler

max_forks_repo_name: gergoerdi/clash-compiler


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE UnboxedTuples #-}

-- |
-- Module      : PSO.Neural
-- Description : Particle Swarm Optimisation
-- Copyright   : (c) Tom Westerhout, 2017
-- License     : BSD3
-- Maintainer  : t.westerhout@student.ru.nl
-- Stability   : experimental
module PSO.Neural
  (
  -- mutableMcmcLoop
  -- , energyHH1DOpen
  -- , energyHH1DOpenC
  -- , energyHH1DOpenMKLC
  -- , energyHH1DOpenZ
  -- , energyVarHH1DOpen
  -- , RBM(..)
    Rbm
  , mkRbm
  , uniformRbm
  , IsRbm(..)
  , fromRbm
  , toRbm
  , unsafeFreezeRbm
  , unsafeThawRbm
  , mcmcHeisenberg1D
  -- , MCMC(..)
  -- , unsafeFreezeMcmc
  -- , Mcmc
  -- , HH1DOpen(..)
  -- , listAll
  , Measurement(..)
  , HasMean(..)
  , HasVar(..)
  ) where

import Prelude hiding (map, mapM, zipWith, zipWithM)
import qualified Prelude as Prelude

import Control.Exception(assert)
import Control.Monad (Monad(..), (>=>), (<=<))
import Control.Monad.Primitive
import Control.Monad.ST

import Data.Bits
import Data.Coerce
import Data.Complex
import Data.Functor.Identity
import Data.NumInstances.Tuple
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as MV
import Data.Vector.Fusion.Stream.Monadic (Stream(..), Step(..), SPEC(..))
import qualified Data.Vector.Fusion.Stream.Monadic as Stream
import qualified Data.Vector.Fusion.Bundle.Monadic as Bundle

import Foreign.C.Types
import Foreign.Storable
import Foreign.ForeignPtr

import GHC.Generics (Generic)

import Lens.Micro
import Lens.Micro.TH

import System.IO.Unsafe

import PSO.Internal.Neural
import PSO.Random
import PSO.Swarm

type family RealOf x

type instance RealOf Float = Float

type instance RealOf Double = Double

type instance RealOf CFloat = CFloat

type instance RealOf CDouble = CDouble

type instance RealOf (Complex a) = a

type family RbmCore a

type instance RbmCore Float = RbmC

-- type instance RbmCore (Complex Double) = RbmZ

newtype Rbm a = R (ForeignPtr (RbmCore a)) deriving (Generic)

newtype MRbm s a = MR (ForeignPtr (RbmCore a)) deriving (Generic)

class Storable a => IsRbm a where
  newRbm :: PrimMonad m => Int -> Int -> m (MRbm (PrimState m) a)
  cloneRbm :: PrimMonad m => Rbm a -> m (MRbm (PrimState m) a)
  getWeights :: PrimMonad m => MRbm (PrimState m) a -> m (V.MVector (PrimState m) a, V.MVector (PrimState m) a)
  getVisible :: PrimMonad m => MRbm (PrimState m) a -> m (V.MVector (PrimState m) a, V.MVector (PrimState m) a)
  getHidden  :: PrimMonad m => MRbm (PrimState m) a -> m (V.MVector (PrimState m) a, V.MVector (PrimState m) a)
  sizeWeights :: Rbm a -> Int
  sizeVisible :: Rbm a -> Int
  sizeHidden  :: Rbm a -> Int
  -- debugPrintRbm :: Rbm a -> IO ()

unsafeThawRbm :: PrimMonad m => Rbm a -> m (MRbm (PrimState m) a)
unsafeThawRbm (R p) = return (MR p)

unsafeFreezeRbm :: PrimMonad m => MRbm (PrimState m) a -> m (Rbm a)
unsafeFreezeRbm (MR p) = return (R p)

-- setWeights :: PrimMonad m => MRbm (PrimState m) a -> V.Vector a -> V.Vector a -> m ()
-- setWeights mrbm wR wI = do
--   (mwR, mwI) <- getWeights mrbm
--   V.copy mwR wR
--   V.copy mwI wI
--
-- setVisible :: PrimMonad m => MRbm (PrimState m) a -> V.Vector a -> V.Vector a -> m ()
-- setVisible mrbm aR aI = do
--   (maR, maI) <- getVisible mrbm
--   V.copy maR aR
--   V.copy maI aI
--
-- setHidden :: PrimMonad m => MRbm (PrimState m) a -> V.Vector a -> V.Vector a -> m ()
-- setHidden mrbm bR bI = do
--   (mbR, mbI) <- getHidden mrbm
--   V.copy mbR bR
--   V.copy mbI bI

fromRbm :: (Monad m, IsRbm a) => Rbm a -> Stream.Stream m a
fromRbm rbm = Bundle.elements $ Bundle.fromVectors (runST $ unsafeThawRbm rbm >>= go)
  where
    go :: IsRbm a => MRbm s a -> ST s [V.Vector a]
    go mrbm = do
      (wR, wI) <- getWeights mrbm
      (aR, aI) <- getVisible mrbm
      (bR, bI) <- getHidden mrbm
      Prelude.mapM V.freeze [wR, wI, aR, aI, bR, bI]

toRbm :: (PrimMonad m, IsRbm a) => Int -> Int -> Stream.Stream m a -> m (Rbm a)
toRbm n m (Stream step t) = do
  mrbm <- newRbm n m
  (wR, wI) <- getWeights mrbm
  (aR, aI) <- getVisible mrbm
  (bR, bI) <- getHidden mrbm
  copy SPEC wR (n * m) 0 t >>=
    copy SPEC wI (n * m) 0 >>=
      copy SPEC aR n 0 >>=
        copy SPEC aI n 0 >>=
          copy SPEC bR m 0 >>=
            copy SPEC bI m 0
  unsafeFreezeRbm mrbm
  where
    copy !_ !v !size !i s
      | i < size = step s >>= \r ->
        case r of
          Yield x s' -> MV.write v i x >> copy SPEC v size (i + 1) s'
          Skip    s' -> copy SPEC v size i s'
          Done       -> error "No, we're not done yet!"
      | otherwise = return s

-- | Creates a new Restricted Boltzmann Machine with given weights.
mkRbm :: IsRbm a
      => V.Vector (Complex a) -- ^ Visible bias \(a^N\).
      -> V.Vector (Complex a) -- ^ Hidden bias \(a^M\).
      -> V.Vector (Complex a) -- ^ Weights \(a^{M\times N}\). *Column-major layout is assumed!*.
      -> Rbm a
mkRbm a@(V.length -> nx1) b@(V.length -> mx1) w@(V.length -> mxn)
  | mxn == nx1 * mx1 = runST $ toRbm nx1 mx1 $ Bundle.elements $
    Bundle.fromVectors [ realPart `V.map` w, imagPart `V.map` w
                       , realPart `V.map` a, imagPart `V.map` a
                       , realPart `V.map` b, imagPart `V.map` b
                       ]
  | otherwise = error "mkRbm: Dimensions mismatch."

-- | Creates an RBM with weights distributed uniformly within given intervals.
uniformRbm ::
     (PrimMonad m, IsRbm a, Ord a, UniformDist m (Complex a))
  => Int -- ^ Number of visible nodes.
  -> Int -- ^ Number of hidden nodes.
  -> (a, a) -- ^ Range for visible bias.
  -> (a, a) -- ^ Range for hidden bias.
  -> (a, a) -- ^ Range for weights.
  -> m (Rbm a)
uniformRbm !n !m !(amin, amax) !(bmin, bmax) !(wmin, wmax)
  | n < 0 || m < 0 = error "uniformRbm: Negative dimensions??"
  | amin > amax || bmin > bmax || wmin > wmax =
    error "uniformRbm: Lower bound is per definition not \
          \greater than the upper bound."
  | otherwise = do
    a <- uniformVector n (amin :+ amin, amax :+ amin)
    b <- uniformVector m (bmin :+ bmin, bmin :+ bmax)
    w <- uniformVector (n * m) (wmin :+ wmin, wmax :+ wmax)
    return $ mkRbm a b w

instance IsRbm Float where
  {-# INLINE newRbm #-}
  newRbm n m = unsafeIOToPrim $ MR <$> _RbmC'construct n m
  {-# INLINE cloneRbm #-}
  cloneRbm (R p) = unsafeIOToPrim $ MR <$> _RbmC'clone p
  {-# INLINE getWeights #-}
  getWeights (MR p) = _RbmC'getWeights p
  {-# INLINE getVisible #-}
  getVisible (MR p) = _RbmC'getVisible p
  {-# INLINE getHidden #-}
  getHidden (MR p) = _RbmC'getHidden p
  {-# INLINE sizeVisible #-}
  sizeVisible (R fp) = unsafePerformIO $
    withForeignPtr fp (peek >=> return . _rbmC'sizeVisible)
  {-# INLINE sizeHidden #-}
  sizeHidden (R fp) = unsafePerformIO $
    withForeignPtr fp (peek >=> return . _rbmC'sizeHidden)
  {-# INLINE sizeWeights #-}
  sizeWeights x = sizeHidden x * sizeVisible x

getDims :: IsRbm a => Rbm a -> (Int, Int)
getDims x = (sizeVisible x, sizeHidden x)
{-# INLINE getDims #-}

zipWithM :: (IsRbm a, IsRbm b, IsRbm c, PrimMonad m)
         => (a -> b -> m c) -> Rbm a -> Rbm b -> m (Rbm c)
zipWithM f x y = assert (getDims x == getDims y) $
  toRbm (sizeVisible x) (sizeHidden x) $ Stream.zipWithM f (fromRbm x) (fromRbm y)

zipWith :: forall a b c. (IsRbm a, IsRbm b, IsRbm c)
        => (a -> b -> c) -> Rbm a -> Rbm b -> Rbm c
zipWith f x y = runST $ zipWithM f' x y
  where f' :: a -> b -> ST s c
        f' a b = return $ f a b

mapM :: (IsRbm a, IsRbm b, PrimMonad m)
     => (a -> m b) -> Rbm a -> m (Rbm b)
mapM f x = toRbm (sizeVisible x) (sizeHidden x) $ Stream.mapM f $ fromRbm x

map :: (IsRbm a, IsRbm b)
     => (a -> b) -> Rbm a -> Rbm b
map f x = runST $ mapM (return . f) x

instance Num (Rbm Float) where
  (+) x y = zipWith (+) x y
  (-) x y = zipWith (-) x y
  (*) x y = zipWith (*) x y
  abs x = map abs x
  signum x = map signum x
  negate x = map negate x
  fromInteger = undefined

instance Scalable Float (Rbm Float) where
  scale λ = map (*λ)

instance (PrimMonad m, DeltaWell m Float Float)
  => DeltaWell m Float (Rbm Float) where
    upDeltaWell κ p x = zipWithM (upDeltaWell κ) p x

{-
instance RBM (Complex Double) where
  newRbm nrVis nrHid = unsafeIOToPrim $ liftM MR (_newRbmZ nrVis nrHid)
  cloneRbm (R p) = unsafeIOToPrim $ liftM MR (_cloneRbmZ p)
  setWeights (MR p) w = unsafeIOToPrim $ _setWeightsZ p (coerce w)
  setVisible (MR p) a = unsafeIOToPrim $ _setVisibleZ p (coerce a)
  setHidden (MR p) b = unsafeIOToPrim $ _setHiddenZ p (coerce b)
  sizeVisible (R p) = _sizeVisibleZ p
  sizeHidden (R p) = _sizeHiddenZ p
  debugPrintRbm (R p) = _printRbmZ p

instance AXPBY (Complex Float) (Complex Float) where
  axpby a (R x) b (MR y) = unsafeIOToPrim $ _caxpbyRbmC (coerce a) x (coerce b) y

instance AXPBY (Complex Double) (Complex Double) where
  axpby a (R x) b (MR y) = unsafeIOToPrim $ _zaxpbyRbmZ (coerce a) x (coerce b) y

instance SCAL (Complex Float) (Complex Float) where
  scal a (MR x) = unsafeIOToPrim $ _cscaleRbmC (coerce a) x

instance SCAL Float (Complex Float) where
  scal a (MR x) = unsafeIOToPrim $ _cscaleRbmC (coerce a :+ 0) x

instance SCAL (Complex Double) (Complex Double) where
  scal a (MR x) = unsafeIOToPrim $ _zscaleRbmZ (coerce a) x

instance SCAL Double (Complex Double) where
  scal a (MR x) = unsafeIOToPrim $ _zscaleRbmZ (coerce a :+ 0) x

instance (Num a, AXPBY a a, SCAL a a) => Num (Rbm a) where
  (+) x y = runST $ cloneRbm x >>= \z ->
    axpby (1 :: a) y (1 :: a) z >> unsafeFreezeRbm z
  (-) x y = runST $ cloneRbm x >>= \z ->
    axpby ((-1) :: a) y (1 :: a) z >> unsafeFreezeRbm z
  (*) x y = undefined
  negate x = runST $ cloneRbm x >>= \z -> scal ((-1) :: a) z >> unsafeFreezeRbm z
  fromInteger = undefined
  abs = undefined
  signum = undefined

instance (SCAL λ a) => Scalable λ (Rbm a) where
  scale k x = runST $ cloneRbm x >>= \z -> scal k z >> unsafeFreezeRbm z

instance MCMC (Complex Float) where
  newMcmc (R rbm) spin = unsafeIOToPrim $
    _newMcmcC rbm (coerce spin) >>= \mcmc -> return (MM mcmc rbm)
  logWF (M p _) = coerce $ _logWFC p
  logQuotient1 (M p _) flip1 = coerce $ _logQuotientWF1C p flip1
  logQuotient2 (M p _) flip1 flip2 = coerce $ _logQuotientWF2C p flip1 flip2
  propose1 m@(M p _) = (^2) . magnitude . exp . logQuotient1 m -- _proposeC1 p flip1
  propose2 (M p _) flip1 flip2 = _propose2C p flip1 flip2
  update1 (MM p _) flip1 = unsafeIOToPrim $ _accept1C p flip1
  update2 (MM p _) flip1 flip2 = unsafeIOToPrim $ _accept2C p flip1 flip2
  debugPrintMcmc (M p _) = _printMcmcC p

instance MCMC (Complex Double) where
  newMcmc (R rbm) spin = unsafeIOToPrim $
    _newMcmcZ rbm (coerce spin) >>= \mcmc -> return (MM mcmc rbm)
  logWF (M p _) = coerce $ _logWFZ p
  logQuotient1 (M p _) flip1 = coerce $ _logQuotientWF1Z p flip1
  logQuotient2 (M p _) flip1 flip2 = coerce $ _logQuotientWF2Z p flip1 flip2
  propose1 m@(M p _) = (^2) . magnitude . exp . logQuotient1 m -- _proposeZ1 p flip1
  propose2 (M p _) flip1 flip2 = _propose2Z p flip1 flip2
  update1 (MM p _) flip1 = unsafeIOToPrim $ _accept1Z p flip1
  update2 (MM p _) flip1 flip2 = unsafeIOToPrim $ _accept2Z p flip1 flip2
  debugPrintMcmc (M p _) = _printMcmcZ p

mutableMcmcLoop ::
     (PrimMonad m)
  => (α -> Int -> m β)
  -> (α -> β -> m Bool)
  -> (α -> β -> m ())
  -> α
  -> Int
  -> Producer α m ()
mutableMcmcLoop propose accept update ψ n = do
  proposal <- lift $ propose ψ n
  isAcceptable <- lift $ accept ψ proposal
  yield ψ
  when isAcceptable (lift (update ψ proposal))
  mutableMcmcLoop propose accept update ψ (n + 1)

class MCMC a => HH1DOpen a where
  locEnergyHH1DOpen :: Mcmc a -> a

instance HH1DOpen (Complex Float) where
  locEnergyHH1DOpen (M p _) = coerce (_locEHH1DOpenC p)

instance HH1DOpen (Complex Double) where
  locEnergyHH1DOpen (M p _) = coerce (_locEHH1DOpenZ p)

unpack2Bools :: (Num α, Bits α) => Int -> α -> V.Vector Bool
unpack2Bools n x = V.fromList $ doUnpack n x []
  where doUnpack 0 _ []   = [False]
        doUnpack 0 _ bits = bits
        doUnpack n x bits = doUnpack (n - 1) (x `shiftR` 1) $
                              (testBit x 0) : bits

listAll :: forall m a.
     ( PrimMonad m
     , HH1DOpen a
     , Num a
     )
  => Rbm a -> m (V.Vector a)
listAll rbm =
  let n = sizeVisible rbm
      toSpin True  = 1
      toSpin False = -1
      σs :: [V.Vector a]
      σs = V.map toSpin <$> unpack2Bools n <$> [(0 :: Int) .. (2 ^ n - 1)]
   in V.fromList <$> (newMcmc rbm >=> unsafeFreezeMcmc >=> return . logWF) `mapM` σs

energyHH1DOpen ::
     forall m a λ.
     ( PrimMonad m
     , Ord (RealOf a)
     , Fractional a
     , HH1DOpen a
     , Randomisable m Bool
     , Randomisable m (RealOf a)
     , UniformDist m Int
     -- , MonadIO m
     -- , a ~ Complex λ
     -- , Fractional λ
     -- , Real λ
     -- , Show λ
     )
  => Rbm a -> Int -> Int -> m a
energyHH1DOpen rbm offset steps =
  let n = sizeVisible rbm
      proposer :: MMcmc (PrimState m) a -> Int -> m Int
      proposer _ i = {- return (i `mod` n) -} uniform (0, n - 1)
      acceptor :: MMcmc (PrimState m) a -> Int -> m Bool
      acceptor mx flip = unsafeFreezeMcmc mx >>= \x ->
                            ((<= propose1 x flip) <$> random)
      updater :: MMcmc (PrimState m) a -> Int -> m ()
      updater = update1
      func :: MMcmc (PrimState m) a -> m a
      func = unsafeFreezeMcmc >=> return . locEnergyHH1DOpen
      mmcmc :: m (MMcmc (PrimState m) a)
      mmcmc = randomSpin n >>= newMcmc rbm
      states = mutableMcmcLoop proposer acceptor updater
   in mmcmc >>= \mmcmc' ->
                  (/ fromIntegral steps)
                    <$> P.sum
                  --     P.fold debugAdd
                  --            (toRational 0, toRational 0)
                  --            (\(x, y) -> fromRational (x / toRational steps)
                  --                         :+ fromRational (y / toRational steps))
                             (states mmcmc' 0 >-> P.drop offset
                                              >-> P.mapM func
                                              >-> P.take steps)

energyHH1DOpenC :: forall m.
     ( PrimMonad m
     , Randomisable m Bool
     , Randomisable m Float
     , UniformDist m Int
     )
  => Rbm (Complex Float) -> Int -> Int -> m (Complex Float)
energyHH1DOpenC rbm offset steps =
  let n = sizeVisible rbm
      mmcmc = randomSpin n >>= newMcmc rbm
      ints :: m (V.Vector Int)
      ints = uniformVector (offset + steps) (0, n - 1)
      floats :: m (V.Vector Float)
      floats = V.replicateM (offset + steps) random
   in ints >>= \ is ->
        floats >>= \ fs ->
          mmcmc >>= \ (MM p (RbmC m)) ->
            unsafeIOToPrim $ (_runMcmcC p offset steps is fs) >>= \ x ->
              touchForeignPtr m >> return (coerce x)

energyHH1DOpenZ :: forall m.
     ( PrimMonad m
     , Randomisable m Bool
     , Randomisable m Double
     , UniformDist m Int
     )
  => Rbm (Complex Double) -> Int -> Int -> m (Complex Double)
energyHH1DOpenZ rbm offset steps =
  let n = sizeVisible rbm
      mmcmc = randomSpin n >>= newMcmc rbm
      ints :: m (V.Vector Int)
      ints = uniformVector (offset + steps) (0, n - 1)
      floats :: m (V.Vector Double)
      floats = V.replicateM (offset + steps) random
   in ints >>= \ is ->
        floats >>= \ fs ->
          mmcmc >>= \ (MM p (RbmZ m)) ->
            unsafeIOToPrim $ (_runMcmcZ p offset steps is fs) >>= \ x ->
              touchForeignPtr m >> return (coerce x)
-}

mcmcHeisenberg1D ::
     PrimMonad m => Rbm Float -> Int -> Int -> m (Measurement Float)
mcmcHeisenberg1D (R rbm) offset steps = unsafeIOToPrim $ _RbmC'heisenberg1D rbm offset steps

instance Eq a => Eq (Measurement a) where
  (==) (Measurement _ x) (Measurement _ y) = x == y

instance Ord a => Ord (Measurement a) where
  -- (<=) (Measurement a x) (Measurement b y) = x <= y
  (<=) (Measurement a x) (Measurement b y) = a <= b

instance Num a => Num (Measurement a) where
  (+) (Measurement a b) (Measurement c d) = Measurement (a + c) (b + d)
  (-) (Measurement a b) (Measurement c d) = Measurement (a - c) (b - d)
  (*) (Measurement a b) (Measurement c d) = Measurement (a * c) (b * d)
  negate (Measurement a b) = Measurement (negate a) (negate b)
  fromInteger n = Measurement (fromInteger n) (fromInteger n)
  abs (Measurement a b) = Measurement (abs a) (abs b)
  signum (Measurement a b) = Measurement (signum a) (signum b)

instance Fractional a => Fractional (Measurement a) where
  (/) (Measurement a b) (Measurement c d) = Measurement (a / c) (b / d)
  recip (Measurement a b) = Measurement (recip a) (recip b)

makeFields ''Measurement

{-

debugAdd :: forall λ. (Real λ, Fractional λ, Show λ)
         => (λ, λ) -> (λ, λ) -> (λ, λ)
debugAdd (x, y) (a, b) = unsafePerformIO $ do
  putStrLn $ show x ++ "\t" ++ show y ++ "\t" ++ show a ++ "\t" ++ show b
  return (x + a, y + b)

energyVarHH1DOpen ::
     forall m a.
     ( PrimMonad m
     , Ord a
     , RealFloat a
     , HH1DOpen (Complex a)
     , Randomisable m Bool
     , Randomisable m a
     , UniformDist m Int
     -- , MonadIO m
     , Show a
     )
  => Rbm (Complex a) -> Int -> Int -> m (Measurement a)
energyVarHH1DOpen rbm offset steps =
  let n = sizeVisible rbm
      proposer :: MMcmc (PrimState m) (Complex a) -> Int -> m Int
      proposer _ i = return (i `mod` n) {- uniform (0, n - 1) -}
      acceptor :: MMcmc (PrimState m) (Complex a) -> Int -> m Bool
      acceptor mx flip = unsafeFreezeMcmc mx >>= \x ->
                            ((<= propose1 x flip) <$> random)
      updater :: MMcmc (PrimState m) (Complex a) -> Int -> m ()
      updater = update1
      func :: MMcmc (PrimState m) (Complex a) -> m (a, a)
      func = unsafeFreezeMcmc >=> return . locEnergyHH1DOpen
                              >=> return . (\ x -> (realPart x, (realPart x)^2))
      mmcmc :: m (MMcmc (PrimState m) (Complex a))
      mmcmc = randomSpin n >>= newMcmc rbm
      states = mutableMcmcLoop proposer acceptor updater
   in mmcmc >>= \mmcmc' ->
                  P.fold debugAdd
                         0
                         (\(h, h2) -> let mean = (h  / fromIntegral steps)
                                          var  = (h2 / fromIntegral steps) - mean^2
                                       in MV mean var)
                         (states mmcmc' 0 >-> P.drop offset
                                          >-> P.mapM func
                                          >-> P.take steps)
-}

max_stars_repo_name: twesterhout/tcm-swarm

max_issues_repo_name: twesterhout/tcm-swarm

max_forks_repo_name: twesterhout/tcm-swarm