text:
module Matrix.StateTest
  (
    testsMatrix
  ) where

import Test.HUnit

import Numeric.LinearAlgebra.HMatrix hiding (corr)
import Matrix.State

mat :: Matrix Double
mat = reshape 3 $ fromList $ take 9 [1..]

vec :: Vector Double
vec = fromList $ take 3 [10..]

test1RE :: Test
test1RE = TestCase (assertEqual "replaceElement - first position"
                    (fromList [1,11,12])
                    (replaceElement 0 1 vec))

test2RE :: Test
test2RE = TestCase (assertEqual "replaceElement - last position"
                    (fromList [10,11,1])
                    (replaceElement 2 1 vec))
test3RE :: Test
test3RE = TestCase (assertEqual "replaceElement - in vector"
                    (fromList [10,2,12])
                    (replaceElement 1 2 vec))

test4RE :: Test
test4RE = TestCase (assertEqual "replaceElement - out of bounds low"
                    (fromList [10,11,12])
                    (replaceElement (-1) 2 vec))

test5RE :: Test
test5RE = TestCase (assertEqual "replaceElement - out of bounds high"
                    (fromList [10,11,12])
                    (replaceElement 4 2 vec))

testsReplaceElement :: Test
testsReplaceElement = TestList [  TestLabel "test1" test1RE
                                , TestLabel "test2" test2RE
                                , TestLabel "test3" test3RE
                                , TestLabel "test4" test4RE
                                , TestLabel "test5" test5RE
                                ]

test1RV :: Test
test1RV = TestCase (assertEqual "replaceVector - first row"
                    (reshape 3 $ fromList [10,11,12,4,5,6,7,8,9])
                    (replaceVector 0 vec mat))

test2RV :: Test
test2RV = TestCase (assertEqual "replaceVector - inner row"
                    (reshape 3 $ fromList [1,2,3,10,11,12,7,8,9])
                    (replaceVector 1 vec mat))

test3RV :: Test
test3RV = TestCase (assertEqual "replaceVector - last row"
                    (reshape 3 $ fromList [1,2,3,4,5,6,10,11,12])
                    (replaceVector 2 vec mat))

test4RV :: Test
test4RV = TestCase (assertEqual "replaceVector - index out of bounds - low"
                    mat
                    (replaceVector (-1) vec mat))

test5RV :: Test
test5RV = TestCase (assertEqual "replaceVector - index out of bounds - high"
                    mat
                    (replaceVector 10 vec mat))

testsReplaceVector :: Test
testsReplaceVector = TestList [ TestLabel "test1" test1RV
                              , TestLabel "test2" test2RV
                              , TestLabel "test3" test3RV
                              , TestLabel "test4" test4RV
                              , TestLabel "test5" test5RV
                              ]

test1RME :: Test
test1RME = TestCase (assertEqual "replaceMatrixElement - first row"
                     (reshape 3 $ fromList [0,2,3,4,5,6,7,8,9])
                     (replaceMatrixElement (0,0) 0 mat))

test2RME :: Test
test2RME = TestCase (assertEqual "replaceMatrixElement - middle row"
                     (reshape 3 $ fromList [1,2,3,0,5,6,7,8,9])
                     (replaceMatrixElement (1,0) 0 mat))

test3RME :: Test
test3RME = TestCase (assertEqual "replaceMatrixElement - end first row"
                     (reshape 3 $ fromList [1,2,0,4,5,6,7,8,9])
                     (replaceMatrixElement (0,2) 0 mat))

test4RME :: Test
test4RME = TestCase (assertEqual "replaceMatrixElement - last row"
                     (reshape 3 $ fromList [1,2,3,4,5,6,7,8,0])
                     (replaceMatrixElement (2,2) 0 mat))

test5RME :: Test
test5RME = TestCase (assertEqual "replaceMatrixElement - out of bounds"
                     mat
                     (replaceMatrixElement (0,4) 0 mat))

test6RME :: Test
test6RME = TestCase (assertEqual "replaceMatrixElement - out of bounds"
                     mat
                     (replaceMatrixElement ((-1),0) 0 mat))

testsReplaceMatrixElement :: Test
testsReplaceMatrixElement = TestList [ TestLabel "test1" test1RME
                                     , TestLabel "test2" test2RME
                                     , TestLabel "test3" test3RME
                                     , TestLabel "test4" test4RME
                                     , TestLabel "test5" test5RME
                                     , TestLabel "test6" test6RME]

testsMatrix :: Test
testsMatrix = TestList [ TestLabel "replaceElementTests" testsReplaceElement
                       , TestLabel "replaceVectorTests" testsReplaceVector
                       , TestLabel "replaceMatrixElementTests" testsReplaceMatrixElement]

max_stars_repo_name: eklinkhammer/neural-algorithms

max_issues_repo_name: eklinkhammer/neural-algorithms

max_forks_repo_name: eklinkhammer/neural-algorithms


text:
{-# OPTIONS_GHC -fno-warn-unused-do-bind #-}

module Scheme48.Parsers where

import Scheme48.Types

import Control.Monad
import Numeric
import Data.Ratio
import Data.Complex
import Data.Array
import Text.ParserCombinators.Parsec hiding (spaces)

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

spaces :: Parser ()
spaces = skipMany1 space

-- Strings --

escapedChars :: Parser Char
escapedChars = do
                _ <- char '\\'
                x <- oneOf "\\\"nrt"
                return $ case x of
                            '\\' -> x
                            '"'  -> x
                            'n'  -> '\n'
                            'r'  -> '\r'
                            't'  -> '\t'
                            _    -> error "impossible"

parseString :: Parser LispVal
parseString = try $ do
                _ <- char '"'
                x <- many (noneOf "\"")
                _ <- char '"'
                return $ String x

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

-- Atoms --

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

--- Numbers ---

parseNumber :: Parser LispVal
parseNumber = try parseComplex
           <|> try parseRatio
           <|> try parseFloat
           <|> try parseDigit1
           <|> try parseDigit2
           <|> try parseHex
           <|> try parseOct
           <|> try parseBin

parseDigit1 :: Parser LispVal
parseDigit1 = try $ liftM (Number . read) (many1 digit)

parseDigit2 :: Parser LispVal
parseDigit2 = try $ string "#d" >> parseDigit1

parseHex :: Parser LispVal
parseHex =  try $ liftM (Number . hex2dig) (string "#x" >> many1 hexDigit)

hex2dig :: (Num a, Eq a) => String -> a
hex2dig = fst . head . readHex

parseOct :: Parser LispVal
parseOct = try $ liftM (Number . oct2dig) (string "#o" >> many1 octDigit)

oct2dig :: (Num a, Eq a) => String -> a
oct2dig = fst . head . readOct

parseBin :: Parser LispVal
parseBin = try $ liftM (Number . bin2dig) (string "#b" >> many1 (oneOf "01"))

bin2dig :: String -> Integer
bin2dig = bin2dig' 0

bin2dig' :: Num a => a -> String -> a
bin2dig' d ""     = d
bin2dig' d (x:xs) = let
                        n = case x of
                                '0' -> 0
                                '1' -> 1
                                _   -> error "Impossible"
                        o = 2 * d + n
                    in bin2dig' o xs

-- Floats --

parseFloat :: Parser LispVal
parseFloat = try $ do
  i <- many1 digit
  _ <- char '.'
  d <- many1 digit
  return $ Float $ fst . head $ readFloat $ i ++ "." ++ d

-- Ratios --

parseRatio :: Parser LispVal
parseRatio = try $ do
  n <- many1 digit
  _ <- char '/'
  d <- many1 digit
  return $ Ratio $ (read n) % (read d)

-- Complex --

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

toDouble :: LispVal -> Double
toDouble (Float f) = f
toDouble (Number n) = fromIntegral n
toDouble e = error $ "wrong type: " ++ show e

-- Booleans --

parseBool :: Parser LispVal
parseBool = try $ char '#' >>
  ((char 't' >> return (Bool True)) <|> (char 'f' >> return (Bool False)))

-- Lists --

parseList :: Parser LispVal
parseList = liftM List $ sepBy parseExpr spaces

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

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

parseLists :: Parser LispVal
parseLists = try $ do
  _ <- char '('
  x <- try parseList <|> parseDottedList
  _ <- char ')'
  return x

-- quasiquote --

parseQuasi :: Parser LispVal
parseQuasi = try $ liftM (\x -> List [Atom "quasiquote", x]) (char '`' >> parseExpr)

parseUnQuote :: Parser LispVal
parseUnQuote = try $ liftM (\x -> List [Atom "unquote", x]) (char ',' >> parseExpr)

-- Vector --

parseVector :: Parser LispVal
parseVector = try $ do
  _ <- string "#("
  v <- sepBy parseExpr spaces
  _ <- char ')'
  return $ Vector $ listArray (0, (length v - 1)) v

-- Core Parser --

parseExpr :: Parser LispVal
parseExpr = parseNumber
        <|> parseBool
        <|> parseCharacter
        <|> parseLists
        <|> parseQuasi
        <|> parseUnQuote
        <|> parseVector
        <|> parseString
        <|> parseAtom
        <|> parseQuoted

parseExprs :: String -> Either ParseError LispVal
parseExprs = parse parseExpr "lisp"

ext: hs


text:
module Operators where

import Data.Complex
import Data.List (foldl1')

import BraKet
import LinAlg

import qualified Numeric.LinearAlgebra as L

-- Creation and annihilation operator definitions --

-- Creation operator in Fock basis
create :: Int -> MatOp
create n = Operator $ (L.><) n n elements
  where
    elements = [if i `elem` nonzeros then (sqrt .fromIntegral) (i `mod` n + 1) else 0 | i <- [0..n^2-1]]
    nonzeros = [i * n + i - 1 | i <- [0..n^2]]

-- Annihilation operator is just conjugate of creation
annihilate :: Int -> MatOp
annihilate n = L.tr <$> create n

-- Number operator definition: n = conj(a) <> a
number :: Int -> MatOp
number n = create n <> annihilate n

-- Density matrix operator
densityMatrix :: [Complex Double] -> [Ket (Complex Double)] -> MatOp
densityMatrix ps kets = foldl1 (<+>) $ zipWith (<**>) ps projectors
    where
      projectors = (\state -> outerProduct state state) <$> kets :: [MatOp]

-- Displacement operator
displacement :: Int -> Complex Double -> MatOp
displacement n alpha = L.expm <$> (conjugate alpha <**> create n
                                        <~> (alpha <**> annihilate n))

squeeze :: Int -> Complex Double -> MatOp
squeeze n z = L.expm <$> ((0.5 :+ 0.0) <**> (conjugate z <**> (annihilate n <> annihilate n)
                                                  <~> (z <**> (create n <> create n))))

max_stars_repo_name: Jvanrhijn/optiqs

max_issues_repo_name: Jvanrhijn/optiqs

max_forks_repo_name: Jvanrhijn/optiqs


text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Trustworthy #-}

{- |
Module      :  Physics.Learn.Schrodinger1D
Copyright   :  (c) Scott N. Walck 2015-2018
License     :  BSD3 (see LICENSE)
Maintainer  :  Scott N. Walck <walck@lvc.edu>
Stability   :  experimental

This module contains functions to
solve the (time dependent) Schrodinger equation
in one spatial dimension for a given potential function.
-}

module Physics.Learn.Schrodinger1D
    (
    -- * Potentials
      freeV
    , harmonicV
    , squareWell
    , doubleWell
    , stepV
    , wall
    -- * Initial wavefunctions
    --    , harm
    , coherent
    , gaussian
    , movingGaussian
    -- * Utilities
    , stateVectorFromWavefunction
    , hamiltonianMatrix
    , expectX
    , picture
    , xRange
    , listForm
    )
    where

import Data.Complex
    ( Complex(..)
    , magnitude
    )
import Graphics.Gloss
    ( Picture(..)
    , yellow
    , black
    , Display(..)
    , display
    )
-- import Math.Polynomial.Hermite
--     ( evalPhysHermite
--     )
import Numeric.LinearAlgebra
    ( R
    , C
    , Vector
    , Matrix
    , (|>)
    , (<.>)
    , fromLists
    , toList
    , size
    )
import Physics.Learn.QuantumMat
    ( probVector
    , timeEv
    )

--i :: Complex Double
--i = 0 :+ 1

----------------
-- Potentials --
----------------

-- | Free potential.
--   The potential energy is zero everywhere.
freeV
    :: Double  -- ^ position
    -> Double  -- ^ potential energy
freeV _x = 0

-- | Harmonic potential.
--   This is the potential energy of a linear spring.
harmonicV
    :: Double  -- ^ spring constant
    -> Double  -- ^ position
    -> Double  -- ^ potential energy
harmonicV k x = k * x**2 / 2

-- | A double well potential.
--   Potential energy is a quartic function of position
--   that gives two wells, each approximately harmonic
--   at the bottom of the well.
doubleWell
    :: Double  -- ^ width (for both wells and well separation)
    -> Double  -- ^ energy height of barrier between wells
    -> Double  -- ^ position
    -> Double  -- ^ potential energy
doubleWell a v0 x = v0 * ((x**2 - a**2)/a**2)**2

-- | Finite square well potential.
--   Potential is zero inside the well,
--   and constant outside the well.
--   Well is centered at the origin.
squareWell
    :: Double  -- ^ well width
    -> Double  -- ^ energy height of well
    -> Double  -- ^ position
    -> Double  -- ^ potential energy
squareWell l v0 x
    | abs x < l/2  = 0
    | otherwise    = v0

-- | A step barrier potential.
--   Potential is zero to left of origin.
stepV
    :: Double  -- ^ energy height of barrier (to the right of origin)
    -> Double  -- ^ position
    -> Double  -- ^ potential energy
stepV v0 x
    | x < 0      = 0
    | otherwise  = v0

-- | A potential barrier with thickness and height.
wall
    :: Double  -- ^ thickness of wall
    -> Double  -- ^ energy height of barrier
    -> Double  -- ^ position of center of barrier
    -> Double  -- ^ position
    -> Double  -- ^ potential energy
wall w v0 x0 x
    | abs (x-x0) < w/2  = v0
    | otherwise         = 0

---------------------------
-- Initial wavefunctions --
---------------------------

-- -- | Harmonic oscillator stationary state
-- harm :: Int          -- ^ nonnegative integer n identifying stationary state
--      -> Double       -- ^ x / sqrt(hbar/(m * omega)), i.e. position
--                      --   in units of sqrt(hbar/(m * omega))
--      -> C            -- ^ complex amplitude
-- harm n u
--     = exp (-u**2/2) * evalPhysHermite n u / sqrt (2^n * fact n * sqrt pi) :+ 0

coherent
    :: R       -- ^ length scale = sqrt(hbar / m omega)
    -> C       -- ^ parameter z
    -> R -> C  -- ^ wavefunction
coherent l z x
    = ((1/(pi*l**2))**0.25 * exp(-x**2/(2*l**2)) :+ 0)
      * exp(-z**2/2 + (sqrt(2/l**2) * x :+ 0) * z)

gaussian
    :: R       -- ^ width parameter
    -> R       -- ^ center of wave packet
    -> R -> C  -- ^ wavefunction
gaussian a x0 x = exp(-(x-x0)**2/(2*a**2)) / sqrt(a * sqrt pi) :+ 0

movingGaussian
    :: R       -- ^ width parameter
    -> R       -- ^ center of wave packet
    -> R       -- ^ l0 = hbar / p0
    -> R -> C  -- ^ wavefunction
movingGaussian a x0 l0 x = exp((0 :+ x/l0) - ((x-x0)**2/(2*a**2) :+ 0)) / (sqrt(a * sqrt pi) :+ 0)

---------------
-- Utilities --
---------------

fact :: Int -> Double
fact 0 = 1
fact n = fromIntegral n * fact (n-1)

linspace :: Double -> Double -> Int -> [Double]
linspace left right num
    = let dx = (right - left) / fromIntegral (num - 1)
      in [ left + dx * fromIntegral n | n <- [0..num-1]]

-- | Transform a wavefunction into a state vector.
stateVectorFromWavefunction :: R         -- ^ lowest x
                            -> R         -- ^ highest x
                            -> Int       -- ^ dimension of state vector
                            -> (R -> C)  -- ^ wavefunction
                            -> Vector C  -- ^ state vector
stateVectorFromWavefunction left right num psi
    = (num |>) [psi x | x <- linspace left right num]

hamiltonianMatrix :: R         -- ^ lowest x
                  -> R         -- ^ highest x
                  -> Int       -- ^ dimension of state vector
                  -> R         -- ^ hbar
                  -> R         -- ^ mass
                  -> (R -> R)  -- ^ potential energy function
                  -> Matrix C  -- ^ Hamiltonian Matrix
hamiltonianMatrix xmin xmax num hbar m pe
    = let coeff = -hbar**2/(2*m)
          dx = (xmax - xmin) / fromIntegral (num - 1)
          diagKEterm = -2 * coeff / dx**2
          offdiagKEterm = coeff / dx**2
          xs = linspace xmin xmax num
      in fromLists [[case abs(i-j) of
                       0  -> (diagKEterm + pe x) :+ 0
                       1  -> offdiagKEterm :+ 0
                       _  -> 0
                          | j <- [1..num] ] | (i,x) <- zip [1..num] xs]

expectX :: Vector C  -- ^ state vector
        -> Vector R  -- ^ vector of x values
        -> R         -- ^ <X>, expectation value of X
expectX psi xs = probVector psi <.> xs

glossScaleX :: Int -> (Double,Double) -> Double -> Float
glossScaleX screenWidth (xmin,xmax) x
    = let w = fromIntegral screenWidth :: Double
      in realToFrac $ (x - xmin) / (xmax - xmin) * w - w / 2

glossScaleY :: Int -> (Double,Double) -> Double -> Float
glossScaleY screenHeight (ymin,ymax) y
    = let h = fromIntegral screenHeight :: Double
      in realToFrac $ (y - ymin) / (ymax - ymin) * h - h / 2

glossScalePoint :: (Int,Int)        -- ^ (screenWidth,screenHeight)
                -> (Double,Double)  -- ^ (xmin,xmax)
                -> (Double,Double)  -- ^ (ymin,ymax)
                -> (Double,Double)  -- ^ (x,y)
                -> (Float,Float)
glossScalePoint (screenWidth,screenHeight) xMinMax yMinMax (x,y)
    = (glossScaleX screenWidth  xMinMax x
      ,glossScaleY screenHeight yMinMax y)

-- | Produce a gloss 'Picture' of state vector
--   for 1D wavefunction.
picture :: (Double, Double)    -- ^ y range
        -> [Double]            -- ^ xs
        -> Vector C            -- ^ state vector
        -> Picture
picture (ymin,ymax) xs psi
    = Color
      yellow
      (Line
       [glossScalePoint
        (screenWidth,screenHeight)
        (head xs, last xs)
        (ymin,ymax)
        p | p <- zip xs (map magSq $ toList psi)])
           where
             magSq = \z -> magnitude z ** 2
             screenWidth  = 1000
             screenHeight =  750

-- options for representing wave functions
-- 1.  A function R -> C
-- 2.  ([R],Vector C), where lengths match
-- 3.  [(R,C)]
-- 4.  (R,R,Vector C)  -- xmin, xmax, state vector (assumes even spacing)

-- 2,4 are best for evolution

listForm :: (R,R,Vector C) -> ([R],Vector C)
listForm (xmin,xmax,v)
    = let dt = (xmax - xmin) / fromIntegral (size v - 1)
      in ([xmin, xmin + dt .. xmax],v)

{-
-- | Given an initial state vector and
--   state propagation function, produce a simulation.
--   The 'Float' in the state propagation function is the time
--   interval for one timestep.
simulate1D :: [Double] -> Vector C -> (Float -> (Float,[Double],Vector C) -> (Float,[Double],Vector C)) -> IO ()
simulate1D xs initial statePropFunc
    = simulate display black 10 (0,initial) displayFunc (const statePropFunc)
      where
        display = InWindow "Animation" (screenWidth,screenHeight) (10,10)
        displayFunc (_t,v) = Color yellow (Line [(

      white (\tFloat -> Pictures [Color blue (Line (points (realToFrac tFloat)))
                                 ,axes (screenWidth,screenHeight) (xmin,xmax) (ymin,ymax)])

-- | Produce a state propagation function from a time-dependent Hamiltonian.
--   The float is dt.
statePropGloss :: (Double -> Matrix C) -> Float -> (Float,Vector C) -> (Float,Vector C)
statePropGloss ham dt (tOld,v)
    = (tNew, timeEv (realToFrac dt) (ham tMid) v)
      where
        tNew = tOld + dt
        tMid = realToFrac $ (tNew + tOld) / 2

-- | Given an initial state vector and a time-dependent Hamiltonian,
--   produce a visualization of a 1D wavefunction.
evolutionBlochSphere :: Vector C -> (Double -> Matrix C) -> IO ()
evolutionBlochSphere psi0 ham
    = simulateBlochSphere 0.01 psi0 (stateProp ham)

-}

{-
def triDiagMatrixMult(square_arr,arr):
    num = len(arr)
    result = array([0 for n in range(num)],dtype=complex128)
    result[0] = square_arr[0][0] * arr[0] + square_arr[0][1] * arr[1]
    for n in range(1,num-1):
        result[n] = square_arr[n][n-1] * arr[n-1] + square_arr[n][n] * arr[n] \
            + square_arr[n][n+1] * arr[n+1]
    result[num-1] = square_arr[num-1][num-2] * arr[num-2] \
        + square_arr[num-1][num-1] * arr[num-1]
    return result
-}

------------------
-- Main program --
------------------

-- n is number of points
-- n-1 is number of intervals
xRange :: R -> R -> Int -> [R]
xRange xmin xmax n
    = let dt = (xmax - xmin) / fromIntegral (n - 1)
      in [xmin, xmin + dt .. xmax]

{-
if __name__ == '__main__':
    m = 1
    omega = 10
    xmin = -2.0
    xmax =  2.0
    num = 256
    num = 128
    dt = 0.0002
    dt = 0.01
    xs = linspace(xmin,xmax,num)
    dx = xs[1] - xs[0]

    super = lambda x: (harm0(m,omega)(x) + harm1(m,omega)(x))/sqrt(2)
    shiftedHarm = lambda x: harm0(m,omega)(x-1)
    coh = coherent(m,omega,1)

    print sum(conj(psi)*psi)*dx

    harmV = harmonicV(m * omega**2)

    V = doubleWell(1,0.1*hbar*omega)
    V = squareWell(1.0,hbar*omega)
    V = harmonicV(m*omega**2)
    V = stepV(10*hbar*omega)
    V = wall(0.1,14.0*hbar*omega,0)
    V = freeV

    H = matrixH(m,xmin,xmax,num,V)
    I = matrixI(num)

    (vals,vecs) = eigh(H)

    E0 = vals[0]
    E1 = vals[1]
    psi0 = normalize(transpose(vecs)[0],dx)
    psi1 = normalize(transpose(vecs)[1],dx)

    psi = func2psi(gaussian(0.3,1),xmin,xmax,num)
    psi = func2psi(coh,xmin,xmax,num)
    psi = func2psi(movingGaussian(0.3,10,-1),xmin,xmax,num)

    psi = psi0
    psi = psi1
    psi = (psi0 + psi1)/sqrt(2)

    E = sum(conj(psi)*triDiagMatrixMult(H,psi)).real*dx

    Escale = hbar*omega

    print E
    print Escale

    leftM  = I + 0.5 * i * H / hbar * dt
    rightM = I - 0.5 * i * H / hbar * dt

    box = display(title='Schrodinger Equation',width=1000,height=1000)

    c = curve(pos = psi2rho(psi,xs))
    c.color = color.blue
    c.radius = 0.02

    ball = sphere(radius=0.05,color=color.red,pos=(expectX(psi,xs),0,0))

    pot_curve = [(x,V(x)/Escale,0) for x in xs if V(x)/Escale < xmax]
    pot = curve(color=color.green,pos=pot_curve,radius=0.01)

    Eline = curve(color=(1,1,0),pos=[(x,E/Escale) for x in xs])
    axis = curve(color=color.white,pos=[(x,0) for x in xs])

    while 1:
        psi = solve(leftM,triDiagMatrixMult(rightM,psi))
        c.pos = psi2rho(psi,xs)
        ball.x = expectX(psi,xs)

To Do:
add combinators for potentials
to shift horizontally and vertically,
and to add potentials

-}

-- Are we committed to SI units for hbar?  No.
-- harmonic oscillator functions depend only on sqrt(hbar/m omega)
-- which is a length parameter
-- for moving gaussian, could give hbar/p0 instead of p0
-- (is that debrogie wavelength?  I think it's h/p0)

max_stars_repo_name: walck/learn-physics

max_issues_repo_name: walck/learn-physics

max_forks_repo_name: walck/learn-physics


text:
{-# LANGUAGE DataKinds #-}

module Matrix
  ( Matrix4
  , Vector3
  , vec3
  , perspective
  , identity
  , withFloatMatrix
  , translate
  , projection
  , normalize
  , xHat
  , yHat
  , zHat
  , cross
  , vector3ToFloatList
  , unwrap
  ) where

import           Data.Maybe
import           Foreign.C.Types
import           Foreign.Ptr
import qualified Numeric.LinearAlgebra        as L
import qualified Numeric.LinearAlgebra.Devel  as D
import           Numeric.LinearAlgebra.Static

type Matrix4 = L 4 4

type Vector3 = R 3

xHat = vec3 1 0 0

yHat = vec3 0 1 0

zHat = vec3 0 0 1

--type Matrix4 = Matrix4 Double
perspective :: Double -> Double -> Double -> Double -> Matrix4
perspective near far fov aspect = build gen
  where
    gen i j =
      case (i, j) of
        (0, 0) -> 1.0 / (aspect * tan (fov / 2.0))
        (1, 1) -> 1.0 / tan (fov / 2.0)
        (2, 2) -> -(far + near) / (far - near)
        (2, 3) -> -2.0 * far * near / (far - near)
        (3, 2) -> -1.0
        _      -> 0.0

identity :: Matrix4
identity = eye

translate :: Vector3 -> Matrix4
translate v = t ||| col (v & 1.0)
  where
    (t, _) = splitCols eye

projection :: Double -> Matrix4
projection = perspective 0.1 10000.0 (pi / 4.0)

withFloatMatrix ::
     Matrix4 -> (D.MatrixOrder -> CInt -> CInt -> Ptr Float -> IO r) -> IO r
withFloatMatrix m f = D.applyRaw singleM id (f (D.orderOf singleM))
  where
    singleM = L.single $ unwrap m

withFloatVector :: Vector3 -> (CInt -> Ptr Float -> IO r) -> IO r
withFloatVector m = D.applyRaw singleV id
  where
    singleV = L.single $ unwrap m

--normalize :: Vector3 -> Vector3
normalize v = realToFrac a * v
  where
    a = 1.0 / norm_2 v

vector3ToFloatList :: Vector3 -> [Float]
vector3ToFloatList = L.toList . L.single . unwrap

max_stars_repo_name: aelg/haskell-render

max_issues_repo_name: aelg/haskell-render

max_forks_repo_name: aelg/haskell-render


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict #-}
module STC.Reversal where

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           STC.DFTArray
import           Types
import FourierPinwheel.Array

{-# INLINE computeReversalR2Z2T0S0 #-}
computeReversalR2Z2T0S0 ::
     (R.Source r (Complex Double))
  => [Double]
  -> R.Array r DIM6 (Complex Double)
  -> R.Array D DIM6 (Complex Double)
computeReversalR2Z2T0S0 theta0Freqs sourceArr =
  R.traverse2
    sourceArr
    (fromListUnboxed (Z :. L.length theta0Freqs) theta0Freqs)
    const $ \f1 f2 idx@(Z :. _ :. _ :. l :. _ :. i :. j) ->
    f1 idx * exp (0 :+ f2 (Z :. l) * pi)

{-# INLINE timeReversal #-}
timeReversal :: DFTArray -> DFTArray
timeReversal arr@(DFTArray _ _ thetaFreqs rFreqs _) =
  repaToDFTArray thetaFreqs rFreqs .
  computeS .
  R.traverse2
    (dftArrayToRepa arr)
    (fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs)
    const $ \fArr fFreq idx@(Z :. _ :. theta :. _ :. _) ->
    fArr idx * (cis $ (fFreq (Z :. theta)) * pi)

{-# INLINE timeReversal' #-}
timeReversal' :: DFTArray -> DFTArray
timeReversal' arr@(DFTArray rows cols thetaFreqs _ _) =
  repaToDFTArray thetaFreqs [] .
  computeS .
  R.traverse2
    (fromUnboxed (Z :. (1 :: Int) :. (L.length thetaFreqs) :. cols :. rows) .
     VS.convert . VS.concat . getDFTArrayVector $
     arr)
    (fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs)
    const $ \fArr fFreq idx@(Z :. _ :. theta :. _ :. _) ->
    fArr idx * (cis $ (fFreq (Z :. theta)) * pi)

{-# INLINE timeReversal1 #-}
timeReversal1 :: Int -> Int -> DFTArray -> DFTArray
timeReversal1 numThetaFreq numRFreq arr@(DFTArray _ _ thetaFreqs rFreqs _) =
  repaToDFTArray thetaFreqs rFreqs .
  computeS .
  R.traverse2
    (dftArrayToRepa1 numThetaFreq numRFreq arr)
    (fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs)
    const $ \fArr fFreq idx@(Z :. _ :. theta :. _ :. _) ->
    fArr idx * (cis $ (fFreq (Z :. theta)) * pi)

{-# INLINE timeReversalRepa #-}
timeReversalRepa ::
     (Num e, RealFloat e, R.Source s (Complex e), Unbox e)
  => [e]
  -> R.Array s DIM4 (Complex e)
  -> R.Array D DIM4 (Complex e)
timeReversalRepa thetaFreqs arr =
  R.traverse2 arr (fromListUnboxed (Z :. L.length thetaFreqs) thetaFreqs) const $ \fArr fFreq idx@(Z :. _ :. theta :. _ :. _) ->
    fArr idx * cis (fFreq (Z :. theta) * pi)

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:
module Main (main) where

import           Debug.Trace

import           Criterion
import           Criterion.IO
import           Criterion.Types
import           Data.Foldable
import           Data.List
import           Data.List.NonEmpty (NonEmpty(..))
import qualified Data.List.NonEmpty as NonEmpty
import           Data.List.Split
import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Statistics.Types
import           System.Environment
import           System.IO
import           Text.Read

main :: IO ()
main = do
  [file] <- getArgs
  res <- readJSONReports file
  case res of
    Left err -> putStrLn err
    Right (_, _, reports) ->
      case traverse process reports of
        Nothing -> putStrLn "Failed to process result"
        Just rs -> do
          let plotData = aggregate rs
          for_ plotData $ \p -> do
            let fileName = plotName p <> ".svg"
            putStrLn ("Generating " <> fileName)
            toFile (def & fo_format .~ SVG) fileName (generatePlot p)

-- We assume that

data Result = Result
  { operation :: String
  , implementation :: String
  , size :: Double
  , iters :: Estimate ConfInt Double
  } deriving Show

process :: Report -> Maybe Result
process r = Result op impl <$> readMaybe sizeStr <*> getIters r
  where
    (sizeStr:impl:rest) = reverse (splitOn "/" (reportName r))
    op = intercalate "/" (reverse rest)

getIters :: Report -> Maybe (Estimate ConfInt Double)
getIters r =
  Map.lookup "iters" . regCoeffs =<<
  find (\r -> regResponder r == "time") (anRegress (reportAnalysis r))

data PlotData = PlotData
  { plotName :: String
  , lines :: [LineData]
  } deriving Show

data LineData = LineData
  { lineName :: String
  , linePoints :: [(Double, Estimate ConfInt Double)]
  } deriving Show

aggregate :: [Result] -> [PlotData]
aggregate rs =
  (map (uncurry PlotData) .
   Map.toList . Map.map (map (uncurry toLineData) . Map.toList))
    byImpl
  where
    byOp :: Map String (NonEmpty Result)
    byOp = aggregateBy operation rs
    byImpl :: Map String (Map String (NonEmpty Result))
    byImpl = Map.map (aggregateBy implementation . toList) byOp
    toLineData :: String -> NonEmpty Result -> LineData
    toLineData k rs = LineData k (map (\r -> (size r, iters r)) (toList rs))

aggregateBy :: Ord k => (a -> k) -> [a] -> Map k (NonEmpty a)
aggregateBy f = Map.fromListWith (<>) . map (\x -> (f x, x :| []))

generatePlot :: PlotData -> EC (Layout Double Double) ()
generatePlot (PlotData n ls) = do
  layout_title .= n
  for_ ls $ \(LineData n' ps) -> do
    let ps' = map (\(x, e) -> (LogValue x, estPoint e)) ps
    plot (line n' [ps'])
    plot (points n' ps')

max_stars_repo_name: cocreature/robinhood

max_issues_repo_name: cocreature/robinhood

max_forks_repo_name: cocreature/robinhood


text:
{-# LANGUAGE DataKinds, GADTs, RankNTypes, TypeOperators #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}
module Algebra.Algorithms.Groebner.ZeroDimSpec where
import Algebra.Algorithms.Groebner
import Algebra.Algorithms.ZeroDim
import Algebra.Internal
import Algebra.Internal
import Algebra.Ring.Ideal
import Algebra.Ring.Polynomial
import Algebra.Ring.Polynomial.Quotient

import           Control.Monad
import           Control.Monad.Random
import           Data.Complex
import           Data.Convertible       (convert)
import qualified Data.Matrix            as M
import           Data.Maybe
import qualified Data.Sized.Builtin     as SV
import qualified Data.Vector            as V
import           Numeric.Field.Fraction (Fraction, (%))
import           Test.Hspec
import           Test.Hspec.QuickCheck
import           Test.QuickCheck        hiding (promote)
import           Utils

asGenListOf :: Gen [a] -> a -> Gen [a]
asGenListOf = const

spec :: Spec
spec = parallel $ do
  describe "solveLinear" $ do
    it "solves data set correctly" $
      forM_ linSet $ \set ->
      solveLinear (M.fromLists $ inputMat set) (V.fromList $ inputVec set)
      `shouldBe` Just (V.fromList $ answer set)
    prop "solves any solvable cases" $
      forAll (resize 10 arbitrary) $ \(MatrixCase ms) ->
      let mat = M.fromLists ms :: M.Matrix (Fraction Integer)
      in rank mat == M.ncols mat ==>
           forAll (vector (length $ head ms)) $ \v ->
             let ans = M.getCol 1 $ mat * M.colVector (V.fromList v)
             in solveLinear mat ans == Just (V.fromList v)
    it "cannot solve unsolvable cases" $ do
      pendingWith "need example"
  describe "univPoly" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    prop "produces elimination ideal's monic generator" $ do
      checkForArity [2..4] prop_univPoly
  describe "radical" $ do
    it "really computes radical" $ do
      pendingWith "We can verify correctness by comparing with singular, but it's not quite smart way..."
{-
      checkForArity [2..4] $ \sdim ->
        forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
        stdReduced (generators $ radical ideal) == calcGroebnerBasis (singIdealFun "radical" ideal)
      -- pendingWith "I couldn't formulate the spec for radical without existentials :-("
-}
  describe "fglm" $ modifyMaxSuccess (const 25) $ modifyMaxSize (const 3) $ do
    prop "computes monomial basis" $
      checkForArity [2..4] $ \sdim ->
        case zeroOrSucc sdim of
          IsZero -> error "impossible"
          IsSucc k ->
            withWitness (lneqSucc k) $ withWitness (lneqZero k) $
            withKnownNat k $
            forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
            let base =
                  reifyQuotient (mapIdeal (changeOrder Lex) ideal) $ \ii ->
                  map quotRepr $ fromJust $ standardMonomials' ii
            in stdReduced (snd $ fglm ideal) == stdReduced base
    prop "computes lex base" $ do
      checkForArity [2..4] $ \sdim ->
        case zeroOrSucc sdim of
          IsZero -> error "impossible"
          IsSucc k -> withKnownNat k $
            withWitness (lneqSucc k) $
            withWitness (lneqZero k) $
            forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
            stdReduced (fst $ fglm ideal)
              == stdReduced (calcGroebnerBasisWith Lex ideal)
    prop "returns lex base in descending order" $
      checkForArity [2..4] $ \sdim ->
      case zeroOrSucc sdim of
        IsZero -> error "impossible"
        IsSucc k -> withKnownNat k $
          case (lneqSucc k, lneqZero k) of
            (Witness, Witness) ->
              forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
              isDescending (map leadingMonomial $ fst $ fglm ideal)
  describe "solve'" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    it "solves equation with admissible error" $ do
      checkForArity [2..4] $ prop_isApproximateZero 1e-5 (solve' 1e-5)
  -- describe "solve''" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
  --   it "solves equation with admissible error" $ do
  --     checkForArity [2..4] $ prop_isApproximateZero 1e-5 (solve'' 1e-10)
  describe "solveViaCompanion" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    it "solves equation with admissible error" $ do
      checkForArity [2..4] $ prop_isApproximateZero 1e-5 (solveViaCompanion 1e-5)
  describe "solveM" $ modifyMaxSuccess (const 50) $ modifyMaxSize (const 4) $ do
    prop "solves equation with admissible error" $ \seed ->
      let gen = mkStdGen seed
      in checkForArity [2..4] $ prop_isApproximateZero 1e-10 (\t -> evalRand (solveM t) gen)

isDescending :: Ord a => [a] -> Bool
isDescending xs = and $ zipWith (>=) xs (drop 1 xs)

prop_isApproximateZero :: KnownNat n
                       => Double
                       -> (forall m. ((0 < m) ~ 'True, KnownNat m) =>
                           Ideal (Polynomial (Fraction Integer) m) -> [Sized m (Complex Double)])
                       -> SNat n -> Property
prop_isApproximateZero err solver sn =
  case zeroOrSucc sn of
    IsSucc k -> withKnownNat k $ forAll (zeroDimOf sn) $ \(ZeroDimIdeal ideal) ->
      (case lneqZero k of
        Witness ->
          let anss = solver ideal
              mul r d = convert r * d
          in all (\as -> all ((<err) . magnitude . substWith mul as) $ generators ideal) anss) :: Bool

prop_univPoly :: KnownNat n => SNat n -> Property
prop_univPoly sdim =
  forAll (zeroDimOf sdim) $ \(ZeroDimIdeal ideal) ->
  let ods = enumOrdinal sdim
  in conjoin $ flip map ods $ \nth ->
  let gen = univPoly nth ideal
  in forAll (unaryPoly sdim nth) $ \f ->
  (f `modPolynomial` [gen] == 0) == (f `isIdealMember` ideal)

rank :: (Ord r, Fractional r) => M.Matrix r -> Int
rank mat =
  let Just (u, _, _, _,_, _) = M.luDecomp' mat
  in V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) $ M.getDiag u

data TestSet = TestSet { inputMat :: [[Fraction Integer]]
                       , inputVec :: [Fraction Integer]
                       , answer   :: [Fraction Integer]
                       } deriving (Show, Eq, Ord)

linSet :: [TestSet]
linSet =
  [TestSet
   [[1 ,0 ,0 ,0 ,0 ]
   ,[0 ,(-2) ,(-2) ,(-2) ,(-2) ]
   ,[0 ,0 ,3 % 2,0 ,(-1) % 2]
   ,[0 ,0 ,0 ,0 ,(-5) % 2]
   ,[0 ,1 ,1 ,1 ,1 ]
   ,[0 ,0 ,(-2) ,1 ,(-1) ]
   ]
   [0 ,(-2) ,19 % 5,14 % 5,1 ,0 ]
   [0 ,(-81) % 25,54 % 25,16 % 5,(-28) % 25]
  ]


text:
{-# LANGUAGE DeriveFunctor, DeriveDataTypeable,DeriveGeneric  #-}
module Statistics.Test.Types (
    Test(..)
  , isSignificant
  , TestResult(..)
  , significant
  , PositionTest(..)
  ) where

import Control.DeepSeq  (NFData(..))
import Data.Data (Typeable, Data)
import GHC.Generics

import Statistics.Types (PValue)

-- | Result of hypothesis testing
data TestResult = Significant    -- ^ Null hypothesis should be rejected
                | NotSignificant -- ^ Data is compatible with hypothesis
                  deriving (Eq,Ord,Show,Typeable,Data,Generic)

instance NFData   TestResult

-- | Result of statistical test.
data Test distr = Test
  { testSignificance :: !(PValue Double)
    -- ^ Probability of getting value of test statistics at least as
    --   extreme as measured.
  , testStatistics   :: !Double
    -- ^ Statistic used for test.
  , testDistribution :: distr
    -- ^ Distribution of test statistics if null hypothesis is correct.
  }
  deriving (Eq,Ord,Show,Typeable,Data,Generic,Functor)

instance (NFData   d) => NFData   (Test d) where
  rnf (Test _ _ a) = rnf a

-- | Check whether test is significant for given p-value.
isSignificant :: PValue Double -> Test d -> TestResult
isSignificant p t
  = significant $ p >= testSignificance t

-- | Test type for test which compare positional (mean,median etc.)
--   information of samples.
data PositionTest
  = SamplesDiffer
    -- ^ Test whether samples differ in position. Null hypothesis is
    --   samples are not different
  | AGreater
    -- ^ Test if first sample (A) is larger than second (B). Null
    --   hypothesis is first sample is not larger than second.
  | BGreater
    -- ^ Test if second sample is larger than first.
  deriving (Eq,Ord,Show,Typeable,Data,Generic)

instance NFData   PositionTest

-- | significant if parameter is 'True', not significant otherwise
significant :: Bool -> TestResult
significant True  = Significant
significant False = NotSignificant

max_stars_repo_name: vmchale/statistics

max_issues_repo_name: vmchale/statistics

max_forks_repo_name: vmchale/statistics


text:
{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies, FlexibleInstances, UndecidableInstances, GADTs, TypeOperators, FlexibleContexts #-}

module Physics.Quantum.Algebra where

import Numeric.Field.Class          (Field)
import Numeric.Algebra.Class        (LeftModule((.*)), Multiplicative((*)), Monoidal(zero))
import Numeric.Additive.Class       (Additive((+)))
import Numeric.Algebra.Unital       (Unital(one))

import Prelude hiding ((+), (*), Complex)

import qualified Data.Map as M
import qualified Data.List as L
import Data.Maybe (catMaybes)

import Data.Complex                 (Complex, conjugate)

import qualified Data.Array.Repa as R
import qualified Data.Array.Repa.Eval as R
import qualified Data.Array.Repa.Repr.Unboxed as R
import qualified Data.Array.Repa.Repr.Vector as R

import Data.Functor.Identity (runIdentity)

import qualified Data.Eigen.SparseMatrix as ES
import qualified Data.Eigen.Matrix as E

-- a vector space m over field k
class (Field k, Additive m) => VectorSpace k m | m -> k where
    (*^) :: k -> m -> m

{-
-- lhs type ‘m’ does not determine rhs type ‘k’
instance (Field k, LeftModule k m) => VectorSpace k m where
    (*^) = (.*)
-}

class VectorSpace k m => LieAlgebra k m | m -> k where
    lieBr :: m -> m -> m

class (VectorSpace k a, VectorSpace k b) => DualVectors k a b | a b -> k where
    dual :: a -> b

class VectorSpace k m => InnerProductSpace k m | m -> k where
    (<.>) :: m -> m -> k

-- Endo-morphism of vector space
-- Data structure for operator may be different for different use cases
-- A Map be suitable for sparse operator, matrix for discreet vector indices, and functions for continuous vector indices.
class VectorSpace k v => SimpleOperator k v o | o -> k, o -> v where
    op :: o -> v -> v

-- Once we have simple operator, we can extract the map information from that
-- and get addition and multiplication operations
data OperatorMap v where
    OperatorMap :: VectorSpace k v => (v -> v) -> OperatorMap v

unOperatorMap :: OperatorMap v -> (v -> v)
unOperatorMap (OperatorMap m) = m

extractOperatorMap :: SimpleOperator k v o => o -> OperatorMap v
extractOperatorMap = OperatorMap . op

instance Additive v => Additive (OperatorMap v) where
    (+) (OperatorMap operator1) (OperatorMap operator2) = OperatorMap $ \vec -> (operator1 vec) + (operator2 vec)

instance  Multiplicative (OperatorMap o) where
    (*) (OperatorMap operator1) (OperatorMap operator2) = OperatorMap $ \vec -> operator1 (operator2 vec)

class (DualVectors k cv v, SimpleOperator k v o) => Operator k v cv o | o -> k where
    (|><|) :: v -> cv -> o

(<|>) :: (InnerProductSpace k v, DualVectors k cv v) => cv -> v -> k
(<|>) covec vec = (dual covec) <.> vec

(|>) :: SimpleOperator k v o => o -> v -> v
(|>) = op

(<|) :: (InnerProductSpace k v, DualVectors k cv v, Operator k v cv o) =>  cv -> o -> (v -> k)
(<|) covector operator  = \vector -> covector <|> (operator |> vector)

class IsComplex a where
    complexConjugate :: a -> a

instance Num a => IsComplex (Complex a) where
    complexConjugate = conjugate

----------------------------------
-- Constructions of vector spaces

-- m is a set of vector indices indices
-- May or may not represent orthogonal vectors, write instance of InnerProductSpace separately
-- kets are represented as maps from indices to coefficients
-- sparse kets have finite number of non zero entries
data SparseKet k m where
    SparseKet :: (Field k, Ord m) => M.Map m k -> SparseKet k m

unSparseKet :: SparseKet k m -> M.Map m k
unSparseKet (SparseKet ketMap) = ketMap

data SparseBra k m where
     SparseBra :: (Field k, Ord m) => M.Map m k -> SparseBra k m

unSparseBra :: SparseBra k m -> M.Map m k
unSparseBra (SparseBra braMap) = braMap

instance (Show k, Show m) => Show (SparseKet k m) where
    show ket@(SparseKet map)  = show map

instance (Show k, Show m) => Show (SparseBra k m) where
    show ket@(SparseBra map)  = show map

instance (Ord m, Field k) => Additive (SparseKet k m) where
    (+) (SparseKet a) (SparseKet b) = SparseKet $ M.unionWith (+) a b

instance (Ord m, Field k) => Additive (SparseBra k m) where
    (+) (SparseBra a) (SparseBra b) = SparseBra $ M.unionWith (+) a b

instance (Ord m, Field k) => VectorSpace k (SparseKet k m) where
    (*^) s (SparseKet v) = SparseKet $ M.map (s *) v

instance (Field k, Ord m) => VectorSpace k (SparseBra k m) where
    (*^) s (SparseBra v) = SparseBra $ M.map (s *) v

instance (Field k, IsComplex k, Ord m) => DualVectors k (SparseKet k m) (SparseBra k m) where
    dual (SparseKet a) = SparseBra $ M.map complexConjugate a

instance (Field k, IsComplex k, Ord m) => DualVectors k (SparseBra k m) (SparseKet k m) where
    dual (SparseBra a) = SparseKet $ M.map complexConjugate a

-- very generic mappings
data InfKet k m = InfKet { unContKet :: (m -> k) }
data InfBra k m = InfBra { unContBra :: (m -> k) }

instance Field k => Additive (InfKet k m) where
    (+) (InfKet a) (InfKet b) = InfKet $ \x -> a x + b x

instance Field k => Additive (InfBra k m) where
    (+) (InfBra a) (InfBra b) = InfBra $ \x -> a x + b x

instance Field k => VectorSpace k (InfKet k m) where
    (*^) s (InfKet v) = InfKet $ \x -> s * (v x)

instance Field k => VectorSpace k (InfBra k m) where
    (*^) s (InfBra v) = InfBra $ \x -> s * (v x)

instance (Field k, IsComplex k) => DualVectors k (InfKet k m) (InfBra k m) where
    dual (InfKet a) = InfBra $ \x -> complexConjugate $ a x

instance (Field k, IsComplex k) => DualVectors k (InfBra k m) (InfKet k m) where
    dual (InfBra a) = InfKet $ \x -> complexConjugate $ a x

-- a vector of coefficients at lattice points
data LatticeKet k m = LatticeKet { unLatticeKet :: R.Array R.D R.DIM1 k }
data LatticeBra k m = LatticeBra { unLatticeBra :: R.Array R.D R.DIM1 k }

instance (Show k, Show m) => Show (LatticeKet k m) where
    -- Delayed representation doesn't have show instance
    show (LatticeKet vec) = show $! compute vec
        where compute :: Show k => R.Array R.D R.DIM1 k -> R.Array R.V R.DIM1 k
              compute = R.computeS

instance (Show k, Show m) => Show (LatticeBra k m) where
    -- Delayed representation doesn't have show instance
    show (LatticeBra vec) = show $! compute vec
        where compute :: Show k => R.Array R.D R.DIM1 k -> R.Array R.V R.DIM1 k
              compute = R.computeS

instance Field k => Additive (LatticeKet k m) where
   (+) (LatticeKet a) (LatticeKet b) = LatticeKet $! R.zipWith (+) a b

instance Field k => Additive (LatticeBra k m) where
   (+) (LatticeBra a) (LatticeBra b) = LatticeBra $! R.zipWith (+) a b

instance Field k => VectorSpace k (LatticeKet k m) where
    (*^) s (LatticeKet v) = LatticeKet $! R.map (s *) v

instance Field k => VectorSpace k (LatticeBra k m) where
    (*^) s (LatticeBra v) = LatticeBra $! R.map (s *) v

instance (Field k, IsComplex k) => DualVectors k (LatticeKet k m) (LatticeBra k m) where
    dual (LatticeKet a) = LatticeBra $! R.map complexConjugate a

instance (Field k, IsComplex k) => DualVectors k (LatticeBra k m) (LatticeKet k m) where
    dual (LatticeBra a) = LatticeKet $! R.map complexConjugate a

-- n x 1 matrix
data SparseMatKet k ck m where
    SparseMatKet :: E.Elem k ck => ES.SparseMatrix k ck -> SparseMatKet k ck m

unSparseMatKet :: SparseMatKet k ck m -> ES.SparseMatrix k ck
unSparseMatKet (SparseMatKet mat) = mat

data SparseMatBra k ck m where
    SparseMatBra :: E.Elem k ck => ES.SparseMatrix k ck -> SparseMatBra k ck m

unSparseMatBra :: SparseMatBra k ck m -> ES.SparseMatrix k ck
unSparseMatBra (SparseMatBra mat) = mat

instance (Show k) => Show (SparseMatKet k ck m) where
    show (SparseMatKet mat) = show mat

instance (Show k) => Show (SparseMatBra k ck m) where
    show (SparseMatBra mat) = show mat

instance Additive (SparseMatKet k ck m) where
    (+) (SparseMatKet a) (SparseMatKet b) = SparseMatKet $! ES.add a b

instance Additive (SparseMatBra k ck m) where
    (+) (SparseMatBra a) (SparseMatBra b) = SparseMatBra $! ES.add a b

instance Field k => VectorSpace k (SparseMatKet k ck m) where
    (*^) s (SparseMatKet v) = SparseMatKet $! ES.scale s v

instance Field k => VectorSpace k (SparseMatBra k ck m) where
    (*^) s (SparseMatBra v) = SparseMatBra $! ES.scale s v

-------------------------------------
-- Constructions of operators

-- Operator formed from outerproduct of two vectors
data OuterOperator k v cv where
    OuterOperator :: (InnerProductSpace k v, DualVectors k cv v) => [(v, cv)] -> OuterOperator k v cv

unOuterOperator :: OuterOperator k v cv -> [(v, cv)]
unOuterOperator (OuterOperator outerProducts) = outerProducts

instance ( InnerProductSpace k v
         , DualVectors k cv v
         )
         => SimpleOperator k v (OuterOperator k v cv) where
    op (OuterOperator outerProducts) vec = L.foldr1 (+) $ map (\(ovec,ocovec) -> (ocovec <|> vec) *^ ovec) $ outerProducts

instance ( InnerProductSpace k v
         , DualVectors k cv v
         )
         => Operator k v cv (OuterOperator k v cv) where
    (|><|) vec covec = OuterOperator [(vec, covec)]

instance Additive (OuterOperator k v cv) where
    (+) (OuterOperator o1) (OuterOperator o2) = OuterOperator $ o1 L.++ o2

-- a sparse operator represented as a map from (i,j) to coefficient of |i><j|
-- Note that (|i><j|) |k> = <j|k> |i> and therefore the operation requires inner product to be defined
data SparseOperator k m = SparseOperator { unSparseOperator :: M.Map (m, m) k }

instance ( Field k, Ord m,
           InnerProductSpace k (SparseKet k m),
           DualVectors k (SparseBra k m) (SparseKet k m)
         )
         => SimpleOperator k (SparseKet k m) (SparseOperator k m) where
    op (SparseOperator opMap) ket = SparseKet $! M.fromListWith (+) $!
        map (\(veci, coveci) -> ( veci
                                , (SparseBra $! M.fromListWith (+) $! catMaybes $! [fmap ((,) coveci) $! M.lookup (veci, coveci) opMap]) <|> ket
                                )
            ) $!
            (M.keys opMap)

instance ( Field k, Ord m,
           InnerProductSpace k (SparseKet k m),
           DualVectors k (SparseBra k m) (SparseKet k m)
         )
         => Operator k (SparseKet k m) (SparseBra k m) (SparseOperator k m) where
    (|><|) (SparseKet vecMap) (SparseBra covecMap) = SparseOperator $! M.fromListWith (+) $!
        catMaybes [(,) (i,j) <$> (fmap (*) (M.lookup i vecMap) <*> (M.lookup j covecMap)) | i <- M.keys vecMap, j <- M.keys covecMap ]

instance ( Field k, Ord m,
           InnerProductSpace k (SparseKet k m),
           DualVectors k (SparseBra k m) (SparseKet k m)
         )
         => Additive (SparseOperator k m) where
    (+) (SparseOperator a) (SparseOperator b) = SparseOperator $! M.unionWith (+) a b

instance ( Field k, Ord m,
           InnerProductSpace k (SparseKet k m),
           DualVectors k (SparseBra k m) (SparseKet k m)
         )
         => Multiplicative (SparseOperator k m) where
     (*) (SparseOperator a) (SparseOperator b) = SparseOperator $! M.fromListWith (+) $! multiplyMaps a b []
         where
             -- chaoose one index i and sum (over j and k) |j><i| |i><k|
             multiplyMaps :: ( Field k, Ord m,
                               InnerProductSpace k (SparseKet k m),
                               DualVectors k (SparseBra k m) (SparseKet k m)
                             )
                             => M.Map (m,m) k -> M.Map (m,m) k -> [((m,m), k)] -> [((m,m), k)]
             multiplyMaps mapA mapB oldResult = case null mapA || null mapB of
                 True -> oldResult
                 False ->
                   -- toList uses list fusion so this should not go through entire map
                   let oneElement = L.take 1 $ M.toList $ mapA
                   in case oneElement of
                          [] -> oldResult
                          elementToUse : _ -> let indexToUse = (snd $! fst elementToUse)

                                                  -- <i|i>
                                                  indexKet = SparseKet $! M.singleton indexToUse one
                                                  indexSquaredNorm =  indexKet <.> indexKet

                                                  -- take all keys from mapA which have this index in snd of key
                                                  (mapAToUse, mapAToPass) = M.partitionWithKey (\key _ -> snd key == indexToUse) mapA
                                                  -- take all keys from mapA which have this index in snd of key
                                                  (mapBToUse, mapBToPass) = M.partitionWithKey (\key _ -> fst key == indexToUse) mapB

                                                  -- we can discard second index from mapAToUse and first index from mapBToUse
                                                  -- we also know that the preserved indices will occur uniquely in index
                                                  -- so we can just take the cross product
                                                  listAToUse = M.toList $! mapAToUse
                                                  listBToUse = M.toList $! mapBToUse
                                                  newResult = [((a,b), va * vb * indexSquaredNorm) | ((a,_), va) <- listAToUse,
                                                                                ((_,b), vb) <- listBToUse
                                                                                ]

                                              in multiplyMaps mapAToPass mapBToPass $! newResult ++ oldResult

instance ( Field k, Ord m,
           InnerProductSpace k (SparseKet k m),
           DualVectors k (SparseBra k m) (SparseKet k m)
         )
         => LeftModule k (SparseOperator k m) where
    (.*) scalar (SparseOperator b) = SparseOperator $! M.map (scalar *) b

-- a more general operator represented as a map from j to Map i_n s_n where operator takes |i_n> to s_n |j>
-- this construct may not suffice when there are uncountable such i_n
data InfOperator k m = InfOperator { unInfOperator :: m -> M.Map m k }

instance (Field k) => SimpleOperator k (InfKet k m) (InfOperator k m) where
    op operator = \(InfKet coeffMap) -> InfKet $ \index -> foldr1 (+) $
        map ( \(invIndex, factor) -> factor * coeffMap invIndex ) $
            M.toList $ unInfOperator operator $ index

instance (Field k, Ord m) => Additive (InfOperator k m) where
    (+) (InfOperator op1) (InfOperator op2) = InfOperator $ \index -> M.unionWith (+) (op1 index) (op2 index)

instance (Field k, Ord m) => Multiplicative (InfOperator k m) where
    (*) (InfOperator op1) (InfOperator op2) = InfOperator $ \index ->
            foldr1 (M.unionWith (+)) $
                map (\(i,s) -> M.map (s *) $ op2 i) $
                    M.toList $ op1 index

-- A matrix to transfom LatticeKet vector
data LatticeOperator k m where
    LatticeOperator :: (Field k, R.Elt k, R.Unbox k, Num k) => R.Array R.D R.DIM2 k -> LatticeOperator k m

unLatticeOperator :: LatticeOperator k m -> R.Array R.D R.DIM2 k
unLatticeOperator (LatticeOperator mat) = mat

-- Num k is required because of sumP in method definition
instance Field k => SimpleOperator k (LatticeKet k m) (LatticeOperator k m) where
    op (LatticeOperator matrix) (LatticeKet vector) = LatticeKet $! R.delay $! runIdentity $! mvMultP matrix vector

instance Additive (LatticeOperator k m) where
    (+) (LatticeOperator op1) (LatticeOperator op2) = LatticeOperator $! R.zipWith (+) op1 op2

-- Num k is required because of sumP in method definition
instance Multiplicative (LatticeOperator k m) where
    (*) (LatticeOperator op1) (LatticeOperator op2) = LatticeOperator $! R.delay $! runIdentity $! mmMultP op1 op2

instance Field k => LeftModule k (LatticeOperator k m) where
    (.*) scalar (LatticeOperator matrix) = LatticeOperator $! R.map (scalar *) matrix

-- Taken from Repa wiki, extended it to Additive and Multiplicative structures
mmMultP :: (Monad m, R.Source r k, Additive k, Multiplicative k, R.Elt k, R.Unbox k, Num k)
        => R.Array r R.DIM2 k
        -> R.Array r R.DIM2 k
        -> m (R.Array R.U R.DIM2 k)
mmMultP a b = R.sumP (R.zipWith (*) aRepl bRepl)
    where
      -- Note: transpose transposes the lowest two dimensions
      t = R.transpose b
      aRepl = R.extend (R.Z R.:.R.All R.:.colsB R.:.R.All) a
      bRepl = R.extend (R.Z R.:.rowsA R.:.R.All R.:.R.All) t

      (R.Z R.:.colsA R.:.rowsA) = R.extent a
      (R.Z R.:.colsB R.:.rowsB) = R.extent b

-- Multiply a martix with a vector
-- based on the above matrix matrix multiplication
mvMultP :: (Monad m, R.Source r k, Additive k, Multiplicative k, R.Elt k, R.Unbox k, Num k)
        => R.Array r R.DIM2 k
        -> R.Array r R.DIM1 k
        -> m (R.Array R.U R.DIM1 k)
mvMultP a b = R.sumP (R.zipWith (*) aRepl bRepl)
    where
      aRepl = a
      bRepl = R.extend (R.Z R.:.rowsA R.:.R.All) b
      (R.Z R.:.rowsA R.:.colsA) = R.extent a

-- A matrix to transfom LatticeKet vector
data SparseMatOperator k ck m where
    SparseMatOperator :: E.Elem k ck => ES.SparseMatrix k ck -> SparseMatOperator k ck m

unSparseMatOperator :: SparseMatOperator k ck m -> ES.SparseMatrix k ck
unSparseMatOperator (SparseMatOperator mat) = mat

-- Num k is required because of sumP in method definition
instance Field k => SimpleOperator k (SparseMatKet k ck m) (SparseMatOperator k ck m) where
    op (SparseMatOperator matrix) (SparseMatKet vector) = SparseMatKet $! ES.mul matrix vector

instance Field k => Additive (SparseMatOperator k ck m) where
    (+) (SparseMatOperator op1) (SparseMatOperator op2) = SparseMatOperator $! ES.add op1 op2

instance Field k => Multiplicative (SparseMatOperator k ck m) where
    (*) (SparseMatOperator op1) (SparseMatOperator op2) = SparseMatOperator $! ES.mul op1 op2

instance Field k => LeftModule k (SparseMatOperator k ck m) where
    (.*) scalar (SparseMatOperator matrix) = SparseMatOperator $! ES.scale scalar matrix

max_stars_repo_name: sarthakbagaria/has-sci

max_issues_repo_name: sarthakbagaria/has-sci

max_forks_repo_name: sarthakbagaria/has-sci


text:
{-# language Rank2Types #-}
{-# language ScopedTypeVariables #-}
{-# language DataKinds #-}
{-# language TypeOperators #-}
{-# language TypeFamilies #-}
{-# language FlexibleInstances #-}
{-# language FlexibleContexts #-}
{-# language MultiParamTypeClasses #-}
{-# language PolyKinds #-}
{-# language GADTs #-}
{-# language ConstraintKinds #-}
{-# language GeneralizedNewtypeDeriving #-}

module Feldspar.Software.Verify.Primitive where

import Feldspar.Sugar
import Feldspar.Representation
import Feldspar.Software.Primitive
import Feldspar.Software.Expression
import Feldspar.Software.Representation hiding (Nil)
import Feldspar.Software.Verify.Command
import Feldspar.Verify.Arithmetic

import Feldspar.Verify.Monad (Verify)
import qualified Feldspar.Verify.FirstOrder as FO
import qualified Feldspar.Verify.Monad as V
import qualified Feldspar.Verify.SMT as SMT
import qualified Feldspar.Verify.Abstract as A

import Data.Struct
import qualified Data.Map.Strict as Map

import qualified Control.Monad.RWS.Strict as S

import qualified SimpleSMT as SMT hiding (not, declareFun)

import qualified Language.Embedded.Expression as Imp
import qualified Language.Embedded.Imperative.CMD as Imp

import qualified Data.Bits as Bits
import qualified Data.Complex as Complex
import Data.Constraint hiding (Sub)
import Data.Int
import Data.Word
import Data.Typeable

import Language.Syntactic

import GHC.Stack

--------------------------------------------------------------------------------
-- *
--------------------------------------------------------------------------------

newtype Symbolic a = Symbolic { unSymbolic :: Rat }
  deriving (Eq, Ord, Show, V.TypedSExpr, V.SMTOrd)

instance V.Fresh (Symbolic a)
  where
    fresh = V.freshSExpr

symbCast :: Symbolic a -> Symbolic b
symbCast = Symbolic . unSymbolic

symbFun :: SymbParam a => String -> [Symbolic a] -> Symbolic a
symbFun name (args :: [Symbolic a]) = V.fromSMT $
  SMT.fun (symbType (undefined :: a) ++ "-" ++ name) (map V.toSMT args)

fromComplexConstant :: (RealFrac a, SymbParam b) => Complex.Complex a -> Symbolic b
fromComplexConstant c = symbFun "complex" [real, imag]
  where
    real = Symbolic $ fromRational $ toRational $ Complex.realPart c
    imag = Symbolic $ fromRational $ toRational $ Complex.imagPart c

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

data SymbFloat
data SymbDouble
data SymbComplexFloat
data SymbComplexDouble

class SymbParam a
  where
    symbType :: a -> String

instance SymbParam SymbFloat  where symbType _ = "float"
instance SymbParam SymbDouble where symbType _ = "double"
instance SymbParam SymbComplexFloat  where symbType _ = "cfloat"
instance SymbParam SymbComplexDouble where symbType _ = "cdouble"

instance SymbParam a => Num (Symbolic a)
  where
    fromInteger = Symbolic . fromInteger
    x + y  = symbFun "+" [x, y]
    x - y  = symbFun "-" [x, y]
    x * y  = symbFun "*" [x, y]
    abs    = smtAbs
    signum = smtSignum

instance SymbParam a => Fractional (Symbolic a) where
  fromRational = Symbolic . fromRational
  x / y = symbFun "/" [x, y]

instance SymbParam a => Floating (Symbolic a) where
  pi      = fromRational (toRational pi)
  exp x   = symbFun "exp" [x]
  log x   = symbFun "log" [x]
  sqrt x  = symbFun "sqrt" [x]
  x ** y  = symbFun "pow" [x, y]
  sin x   = symbFun "sin" [x]
  cos x   = symbFun "cos" [x]
  tan x   = symbFun "tan" [x]
  asin x  = symbFun "asin" [x]
  acos x  = symbFun "acos" [x]
  atan x  = symbFun "atan" [x]
  sinh x  = symbFun "sinh" [x]
  cosh x  = symbFun "cosh" [x]
  tanh x  = symbFun "tanh" [x]
  asinh x = symbFun "asinh" [x]
  acosh x = symbFun "acosh" [x]
  atanh x = symbFun "atanh" [x]

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

class Floating a => Complex a
  where
    type RealPart a
    complex   :: RealPart a -> RealPart a -> a
    polar     :: RealPart a -> RealPart a -> a
    real      :: a -> RealPart a
    imag      :: a -> RealPart a
    magnitude :: a -> RealPart a
    phase     :: a -> RealPart a
    conjugate :: a -> a

instance Complex (V.SMTExpr Prim (Complex.Complex Float))
  where
    type RealPart (V.SMTExpr Prim (Complex.Complex Float)) = V.SMTExpr Prim Float
    complex   (Float x) (Float y) = ComplexFloat (complex x y)
    polar     (Float x) (Float y) = ComplexFloat (polar x y)
    real      (ComplexFloat x)    = Float (real x)
    imag      (ComplexFloat x)    = Float (imag x)
    magnitude (ComplexFloat x)    = Float (magnitude x)
    phase     (ComplexFloat x)    = Float (phase x)
    conjugate (ComplexFloat x)    = ComplexFloat (conjugate x)

instance Complex (V.SMTExpr Prim (Complex.Complex Double))
  where
    type RealPart (V.SMTExpr Prim (Complex.Complex Double)) = V.SMTExpr Prim Double
    complex   (Double x) (Double y) = ComplexDouble (complex x y)
    polar     (Double x) (Double y) = ComplexDouble (polar x y)
    real      (ComplexDouble x)     = Double (real x)
    imag      (ComplexDouble x)     = Double (imag x)
    magnitude (ComplexDouble x)     = Double (magnitude x)
    phase     (ComplexDouble x)     = Double (phase x)
    conjugate (ComplexDouble x)     = ComplexDouble (conjugate x)

witnessComplex :: (SoftwarePrimType a, SoftwarePrimType (Complex.Complex a)) =>
  Prim (Complex.Complex a) ->
  Dict ( Floating (V.SMTExpr Prim a)
       , Complex  (V.SMTExpr Prim (Complex.Complex a))
       , RealPart (V.SMTExpr Prim (Complex.Complex a)) ~ V.SMTExpr Prim a)
witnessComplex (_ :: Prim (Complex.Complex a)) =
  case softwareRep :: SoftwarePrimTypeRep (Complex.Complex a) of
    ComplexFloatST  -> Dict
    ComplexDoubleST -> Dict

witnessFractional :: (SoftwarePrimType a, Fractional a) =>
  Prim a ->
  Dict (Floating (V.SMTExpr Prim a))
witnessFractional (_ :: Prim a) = case softwareRep :: SoftwarePrimTypeRep a of
    FloatST         -> Dict
    DoubleST        -> Dict
    ComplexFloatST  -> Dict
    ComplexDoubleST -> Dict

witnessIntegral :: (SoftwarePrimType a, Integral a) =>
  Prim a ->
  Dict (Integral (V.SMTExpr Prim a), Bits (V.SMTExpr Prim a))
witnessIntegral (_ :: Prim a) = case softwareRep :: SoftwarePrimTypeRep a of
  Int8ST   -> Dict
  Int16ST  -> Dict
  Int32ST  -> Dict
  Int64ST  -> Dict
  Word8ST  -> Dict
  Word16ST -> Dict
  Word32ST -> Dict
  Word64ST -> Dict

witnessBits :: (SoftwarePrimType a, Bits.Bits a) =>
  Prim a ->
  Dict ( Num a
       , Integral (V.SMTExpr Prim a)
       , Bits (V.SMTExpr Prim a))
witnessBits (_ :: Prim a) = case softwareRep :: SoftwarePrimTypeRep a of
  Int8ST   -> Dict
  Int16ST  -> Dict
  Int32ST  -> Dict
  Int64ST  -> Dict
  Word8ST  -> Dict
  Word16ST -> Dict
  Word32ST -> Dict
  Word64ST -> Dict

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

toRat :: (SoftwarePrimType a, Fractional a) => V.SMTExpr Prim a -> Rat
toRat (x :: V.SMTExpr Prim a) = case softwareRep :: SoftwarePrimTypeRep a of
  FloatST         -> let Float         (Symbolic y) = x in y
  DoubleST        -> let Double        (Symbolic y) = x in y
  ComplexFloatST  -> let ComplexFloat  (Symbolic y) = x in y
  ComplexDoubleST -> let ComplexDouble (Symbolic y) = x in y

fromRat :: forall a. (SoftwarePrimType a, Num a) => Rat -> V.SMTExpr Prim a
fromRat x = case softwareRep :: SoftwarePrimTypeRep a of Int8ST -> Int8 (f2i x)
  where
    f2i :: forall s w. (Sign s, Width w) => Rat -> BV s w
    f2i (Rat x) = BV $ SMT.List [SMT.fam "int2bv" [width (undefined :: w)], SMT.fun "to_int" [x]]

i2n :: forall a b.
  ( V.SMTEval Prim a, SoftwarePrimType a, Integral a
  , V.SMTEval Prim b, SoftwarePrimType b, Num b) =>
  V.SMTExpr Prim a ->
  V.SMTExpr Prim b
i2n x = toBV x $ case softwareRep :: SoftwarePrimTypeRep b of
  Int8ST   -> Int8   . i2i
  Int16ST  -> Int16  . i2i
  Int32ST  -> Int32  . i2i
  Int64ST  -> Int64  . i2i
  Word8ST  -> Word8  . i2i
  Word16ST -> Word16 . i2i
  Word32ST -> Word32 . i2i
  Word64ST -> Word64 . i2i
  FloatST  -> Float  . Symbolic . i2f
  DoubleST -> Double . Symbolic . i2f
  ComplexFloatST  -> ComplexFloat  . Symbolic . i2f
  ComplexDoubleST -> ComplexDouble . Symbolic . i2f
  where
    toBV :: forall a b. (V.SMTEval Prim a, SoftwarePrimType a, Integral a) =>
      V.SMTExpr Prim a -> (forall s w. (Sign s, Width w) => BV s w -> b) -> b
    toBV (x :: V.SMTExpr Prim a) k = case softwareRep :: SoftwarePrimTypeRep a of
      Int8ST   -> let Int8   y = x in k y
      Int16ST  -> let Int16  y = x in k y
      Int32ST  -> let Int32  y = x in k y
      Int64ST  -> let Int64  y = x in k y
      Word8ST  -> let Word8  y = x in k y
      Word16ST -> let Word16 y = x in k y
      Word32ST -> let Word32 y = x in k y
      Word64ST -> let Word64 y = x in k y

    i2f :: (Sign s, Width w) => BV s w -> Rat
    i2f (BV x) = Rat (SMT.fun "to_real" [SMT.fun "bv2int" [x]])

    i2i :: forall s1 w1 s2 w2. (Sign s1, Width w1, Sign s2, Width w2) => BV s1 w1 -> BV s2 w2
    i2i x = case compare m n of
      LT | isSigned x -> V.fromSMT (SMT.signExtend (n-m) (V.toSMT x))
         | otherwise  -> V.fromSMT (SMT.zeroExtend (n-m) (V.toSMT x))
      EQ -> V.fromSMT (V.toSMT x)
      GT -> V.fromSMT (SMT.extract (V.toSMT x) (n-1) 0)
      where
        m = width (undefined :: w1)
        n = width (undefined :: w2)

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

class SymbParam a => SymbComplex a
  where
    type SymbRealPart a

instance SymbComplex SymbComplexFloat
  where
    type SymbRealPart SymbComplexFloat  = SymbFloat

instance SymbComplex SymbComplexDouble
  where
    type SymbRealPart SymbComplexDouble = SymbDouble

instance SymbComplex a => Complex (Symbolic a)
  where
    type RealPart (Symbolic a) = Symbolic (SymbRealPart a)
    complex x y = symbFun "complex" [symbCast x, symbCast y]
    polar x y   = symbFun "polar" [symbCast x, symbCast y]
    real x      = symbCast (symbFun "real" [x])
    imag x      = symbCast (symbFun "imag" [x])
    magnitude x = symbCast (symbFun "magnitude" [x])
    phase x     = symbCast (symbFun "phase" [x])
    conjugate x = symbFun "conjugate" [x]

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

declareSymbFun :: SymbParam a => String -> a ->
  [SMT.SExpr] -> SMT.SExpr -> SMT.SMT ()
declareSymbFun name (_ :: a) args res =
  S.void $ SMT.declareFun (symbType (undefined :: a) ++ "-" ++ name) args res

declareSymbArith :: SymbParam a => a -> SMT.SMT ()
declareSymbArith x = do
  declareSymbFun "+" x     [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "-" x     [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "*" x     [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "/" x     [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "exp" x   [SMT.tReal] SMT.tReal
  declareSymbFun "log" x   [SMT.tReal] SMT.tReal
  declareSymbFun "sqrt" x  [SMT.tReal] SMT.tReal
  declareSymbFun "pow" x   [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "sin" x   [SMT.tReal] SMT.tReal
  declareSymbFun "cos" x   [SMT.tReal] SMT.tReal
  declareSymbFun "tan" x   [SMT.tReal] SMT.tReal
  declareSymbFun "asin" x  [SMT.tReal] SMT.tReal
  declareSymbFun "acos" x  [SMT.tReal] SMT.tReal
  declareSymbFun "atan" x  [SMT.tReal] SMT.tReal
  declareSymbFun "sinh" x  [SMT.tReal] SMT.tReal
  declareSymbFun "cosh" x  [SMT.tReal] SMT.tReal
  declareSymbFun "tanh" x  [SMT.tReal] SMT.tReal
  declareSymbFun "asinh" x [SMT.tReal] SMT.tReal
  declareSymbFun "acosh" x [SMT.tReal] SMT.tReal
  declareSymbFun "atanh" x [SMT.tReal] SMT.tReal

declareSymbComplex :: SymbParam a => a -> SMT.SMT ()
declareSymbComplex x = do
  declareSymbArith x
  declareSymbFun "complex" x   [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "polar" x     [SMT.tReal, SMT.tReal] SMT.tReal
  declareSymbFun "real" x      [SMT.tReal] SMT.tReal
  declareSymbFun "imag" x      [SMT.tReal] SMT.tReal
  declareSymbFun "magnitude" x [SMT.tReal] SMT.tReal
  declareSymbFun "phase" x     [SMT.tReal] SMT.tReal
  declareSymbFun "conjugate" x [SMT.tReal] SMT.tReal

declareFeldsparGlobals :: SMT.SMT ()
declareFeldsparGlobals = do
  declareSymbArith   (undefined :: SymbFloat)
  declareSymbArith   (undefined :: SymbDouble)
  declareSymbComplex (undefined :: SymbComplexFloat)
  declareSymbComplex (undefined :: SymbComplexDouble)
  SMT.declareFun "skolem-int8"   [] (SMT.tBits 8)
  SMT.declareFun "skolem-int16"  [] (SMT.tBits 16)
  SMT.declareFun "skolem-int32"  [] (SMT.tBits 32)
  SMT.declareFun "skolem-int64"  [] (SMT.tBits 64)
  SMT.declareFun "skolem-word8"  [] (SMT.tBits 8)
  SMT.declareFun "skolem-word16" [] (SMT.tBits 16)
  SMT.declareFun "skolem-word32" [] (SMT.tBits 32)
  SMT.declareFun "skolem-word64" [] (SMT.tBits 64)
  return ()

instance FO.Substitute Prim
  where
    type SubstPred Prim = SoftwarePrimType
    subst sub (Prim exp) = Prim (everywhereUp f exp)
      where
        f :: ASTF SoftwarePrimDomain a -> ASTF SoftwarePrimDomain a
        f (Sym (FreeVar x :&: ty)) = case FO.lookupSubst sub (Imp.ValComp x) of
          Imp.ValComp y -> Sym (FreeVar y :&: ty)
          Imp.ValRun  z -> Sym (Lit z :&: ty)
        f (Sym (ArrIx iarr :&: ty) :$ i) = Sym (ArrIx iarr' :&: ty) :$ i
          where iarr' = FO.lookupSubst sub iarr
        f x = x

instance FO.TypeablePred SoftwarePrimType
  where
    witnessTypeable Dict = Dict

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

instance V.SMTEval Prim Bool where
  fromConstant = Bool . SMT.bool
  witnessOrd _ = Dict

instance V.SMTEval Prim Int8 where
  fromConstant = Int8 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-int8" [])

instance V.SMTEval Prim Int16 where
  fromConstant = Int16 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-int16" [])

instance V.SMTEval Prim Int32 where
  fromConstant = Int32 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-int32" [])

instance V.SMTEval Prim Int64 where
  fromConstant = Int64 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-int64" [])

instance V.SMTEval Prim Word8 where
  fromConstant = Word8 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-word8" [])

instance V.SMTEval Prim Word16 where
  fromConstant = Word16 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-word16" [])

instance V.SMTEval Prim Word32 where
  fromConstant = Word32 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-word32" [])

instance V.SMTEval Prim Word64 where
  fromConstant = Word64 . fromIntegral
  witnessNum _ = Dict
  witnessOrd _ = Dict
  skolemIndex  = V.fromSMT (SMT.fun "skolem-word64" [])

instance V.SMTEval Prim Float where
  fromConstant = Float . fromRational . toRational
  witnessOrd _ = Dict
  witnessNum _ = Dict

instance V.SMTEval Prim Double where
  fromConstant = Double . fromRational . toRational
  witnessOrd _ = Dict
  witnessNum _ = Dict

instance V.SMTEval Prim (Complex.Complex Float) where
  fromConstant = ComplexFloat . fromComplexConstant
  witnessNum _ = Dict

instance V.SMTEval Prim (Complex.Complex Double) where
  fromConstant = ComplexDouble . fromComplexConstant
  witnessNum _ = Dict

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

instance V.SMTEval1 Prim where
  type Pred Prim = SoftwarePrimType
  newtype SMTExpr Prim Bool   = Bool SMT.SExpr
    deriving (Typeable)
  newtype SMTExpr Prim Float  = Float  (Symbolic SymbFloat)
    deriving (Typeable, Num, Fractional, Floating, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Double = Double (Symbolic SymbDouble)
    deriving (Typeable, Num, Fractional, Floating, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim (Complex.Complex Float) =
      ComplexFloat (Symbolic SymbComplexFloat)
    deriving (Typeable, Num, Fractional, Floating, V.TypedSExpr)
  newtype SMTExpr Prim (Complex.Complex Double) =
      ComplexDouble (Symbolic SymbComplexDouble)
    deriving (Typeable, Num, Fractional, Floating, V.TypedSExpr)
  newtype SMTExpr Prim Int8   = Int8   (BV Signed W8)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Int16  = Int16  (BV Signed W16)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Int32  = Int32  (BV Signed W32)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Int64  = Int64  (BV Signed W64)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Word8  = Word8  (BV Unsigned W8)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Word16 = Word16 (BV Unsigned W16)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Word32 = Word32 (BV Unsigned W32)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)
  newtype SMTExpr Prim Word64 = Word64 (BV Unsigned W64)
    deriving (Typeable, Num, Real, Enum, Integral, Bits, V.SMTOrd, V.TypedSExpr)

  eval (Prim exp :: Prim a) =
    simpleMatch (\(exp :&: ty) -> case softwarePrimWitType ty of
      Dict -> case V.witnessPred (undefined :: Prim a) of
        Dict -> verifyPrim exp) exp

  witnessPred (_ :: Prim a) = case softwareRep :: SoftwarePrimTypeRep a of
    BoolST   -> Dict
    Int8ST   -> Dict
    Int16ST  -> Dict
    Int32ST  -> Dict
    Int64ST  -> Dict
    Word8ST  -> Dict
    Word16ST -> Dict
    Word32ST -> Dict
    Word64ST -> Dict
    FloatST  -> Dict
    DoubleST -> Dict
    ComplexFloatST  -> Dict
    ComplexDoubleST -> Dict

instance V.SMTOrd (V.SMTExpr Prim Bool) where
  Bool x .<.  Bool y = SMT.not x SMT..&&. y
  Bool x .>.  Bool y = x SMT..&&. SMT.not y
  Bool x .<=. Bool y = SMT.not x SMT..||. y
  Bool x .>=. Bool y = x SMT..||. SMT.not y

instance V.TypedSExpr (V.SMTExpr Prim Bool) where
  smtType _ = SMT.tBool
  toSMT (Bool x) = x
  fromSMT x = Bool x

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

verifyPrim :: forall a .
  (SoftwarePrimType (DenResult a), V.SMTEval Prim (DenResult a), HasCallStack) =>
  SoftwarePrim a ->
  Args (AST SoftwarePrimDomain) a ->
  V.Verify (V.SMTExpr Prim (DenResult a))
verifyPrim (FreeVar x) _ = peekValue x
verifyPrim (Lit x) _     = return (V.fromConstant x)
verifyPrim Add (x :* y :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    S.liftM2 (+) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Sub (x :* y :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    S.liftM2 (-) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Mul (x :* y :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    S.liftM2 (*) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Neg (x :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    fmap negate (V.eval (Prim x))
verifyPrim Abs (x :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    fmap abs (V.eval (Prim x))
verifyPrim Sign (x :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) =
    fmap signum (V.eval (Prim x))
verifyPrim Div (x :* y :* Nil)
  | Dict <- witnessIntegral (undefined :: Prim (DenResult a)) =
    S.liftM2 div (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Mod (x :* y :* Nil)
  | Dict <- witnessIntegral (undefined :: Prim (DenResult a)) =
    S.liftM2 mod (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Quot (x :* y :* Nil)
  | Dict <- witnessIntegral (undefined :: Prim (DenResult a)) =
    S.liftM2 quot (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Rem (x :* y :* Nil)
  | Dict <- witnessIntegral (undefined :: Prim (DenResult a)) =
    S.liftM2 rem (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim FDiv (x :* y :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    S.liftM2 (/) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Pi Nil
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    return pi
verifyPrim Exp (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap exp (V.eval (Prim x))
verifyPrim Log (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap log (V.eval (Prim x))
verifyPrim Sqrt (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap sqrt (V.eval (Prim x))
verifyPrim Pow (x :* y :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    S.liftM2 (**) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Sin (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap sin (V.eval (Prim x))
verifyPrim Cos (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap cos (V.eval (Prim x))
verifyPrim Tan (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap tan (V.eval (Prim x))
verifyPrim Asin (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap asin (V.eval (Prim x))
verifyPrim Acos (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap acos (V.eval (Prim x))
verifyPrim Atan (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap atan (V.eval (Prim x))
verifyPrim Sinh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap sinh (V.eval (Prim x))
verifyPrim Cosh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap cosh (V.eval (Prim x))
verifyPrim Tanh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap tanh (V.eval (Prim x))
verifyPrim Asinh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap asinh (V.eval (Prim x))
verifyPrim Acosh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap acosh (V.eval (Prim x))
verifyPrim Atanh (x :* Nil)
  | Dict <- witnessFractional (undefined :: Prim (DenResult a)) =
    fmap atanh (V.eval (Prim x))
verifyPrim Complex (x :* y :* Nil)
  | Dict <- witnessComplex (undefined :: Prim (DenResult a)) =
    S.liftM2 complex (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Polar (x :* y :* Nil)
  | Dict <- witnessComplex (undefined :: Prim (DenResult a)) =
    S.liftM2 polar (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim Real ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessComplex (undefined :: Prim b) =
    fmap real (V.eval (Prim x))
verifyPrim Imag ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessComplex (undefined :: Prim b) =
    fmap imag (V.eval (Prim x))
verifyPrim Magnitude ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessComplex (undefined :: Prim b) =
    fmap magnitude (V.eval (Prim x))
verifyPrim Phase ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessComplex (undefined :: Prim b) =
    fmap phase (V.eval (Prim x))
verifyPrim Conjugate ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessComplex (undefined :: Prim b) =
    fmap conjugate (V.eval (Prim x))
verifyPrim I2N ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessNum (undefined :: Prim b) = do
    fmap i2n (V.eval (Prim x))
verifyPrim I2B ((x :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessNum (undefined :: Prim b) = do
    x <- V.eval (Prim x)
    return (V.fromSMT (SMT.not (x V..==. 0)))
verifyPrim B2I (x :* Nil)
  | Dict <- V.witnessNum (undefined :: Prim (DenResult a)) = do
    x <- V.eval (Prim x)
    return (V.smtIte (V.toSMT x) 1 0)
verifyPrim Round ((x :: ASTF SoftwarePrimDomain b) :* Nil) = do
  x <- V.eval (Prim x)
  return (fromRat (toRat x))
verifyPrim Not (x :* Nil) =
  fmap (V.fromSMT . SMT.not . V.toSMT) (V.eval (Prim x))
verifyPrim And (x :* y :* Nil) = do
  x <- V.eval (Prim x)
  y <- V.eval (Prim y)
  return (V.fromSMT (V.toSMT x SMT..&&. V.toSMT y))
verifyPrim Or (x :* y :* Nil) = do
  x <- V.eval (Prim x)
  y <- V.eval (Prim y)
  return (V.fromSMT (V.toSMT x SMT..||. V.toSMT y))
verifyPrim Eq ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (V..==.) (V.eval (Prim x)) (V.eval (Prim y)))
verifyPrim Neq ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (./=.) (V.eval (Prim x)) (V.eval (Prim y)))
  where
    x ./=. y = SMT.not (x V..==. y)
verifyPrim Lt ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessOrd  (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (V..<.) (V.eval (Prim x)) (V.eval (Prim y)))
verifyPrim Gt ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessOrd  (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (V..>.) (V.eval (Prim x)) (V.eval (Prim y)))
verifyPrim Lte ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessOrd  (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (V..<=.) (V.eval (Prim x)) (V.eval (Prim y)))
verifyPrim Gte ((x :: ASTF SoftwarePrimDomain b) :* y :* Nil)
  | Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- V.witnessOrd  (undefined :: Prim b) =
    fmap V.fromSMT (S.liftM2 (V..>=.) (V.eval (Prim x)) (V.eval (Prim y)))
verifyPrim BitAnd (x :* y :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)) =
    S.liftM2 (.&.) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim BitOr (x :* y :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)) =
    S.liftM2 (.|.) (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim BitXor (x :* y :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)) =
    S.liftM2 xor (V.eval (Prim x)) (V.eval (Prim y))
verifyPrim BitCompl (x :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)) =
    fmap complement (V.eval (Prim x))
verifyPrim ShiftL (x :* (y :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)),
    Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- witnessIntegral (undefined :: Prim b) = do
    -- todo: should check for undefined behaviour
    x <- V.eval (Prim x)
    y <- V.eval (Prim y)
    return (shiftL x (i2n y))
verifyPrim ShiftR (x :* (y :: ASTF SoftwarePrimDomain b) :* Nil)
  | Dict <- witnessBits (undefined :: Prim (DenResult a)),
    Dict <- V.witnessPred (undefined :: Prim b),
    Dict <- witnessIntegral (undefined :: Prim b) = do
    -- todo: should check for undefined behaviour
    x <- V.eval (Prim x)
    y <- V.eval (Prim y)
    return (shiftR x (i2n y))
verifyPrim (ArrIx (Imp.IArrComp name :: Imp.IArr Index b)) (i :* Nil) = do
  i <- V.eval (Prim i)
  readArray name i
verifyPrim Cond (cond :* x :* y :* Nil) =
  S.liftM3 V.smtIte
    (fmap V.toSMT (V.eval (Prim cond)))
    (V.eval (Prim x))
    (V.eval (Prim y))
verifyPrim exp _ = error ("Unimplemented: " ++ show exp)

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

max_stars_repo_name: markus-git/co-feldspar

max_issues_repo_name: markus-git/co-feldspar

max_forks_repo_name: markus-git/co-feldspar


text:
{-# LANGUAGE FlexibleContexts, FlexibleInstances #-}
module Test.Utils where
import Numeric.LinearAlgebra
import Control.Exception
import Control.Monad
import qualified Data.NeuralNetwork.Backend.BLASHS.Utils as U
import qualified Numeric.LinearAlgebra as L
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector as BV
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as MV
import System.IO.Unsafe

asHM (U.DenseMatrix r c v) = L.reshape c $ unsafePerformIO $ V.freeze v
asDM m = let (r,c) = size m in U.DenseMatrix r c (unsafePerformIO $ V.thaw $ L.flatten m)

good_corr2 :: Int -> L.Matrix Float -> L.Matrix Float -> L.Matrix Float

good_corr2 p k m | w > s     = good_corr2 p m k
                 | otherwise = corr2 k padded
  where
    (w,h) = L.size k
    (s,t) = L.size m
    padded = fromBlocks [[z,0,0]
                        ,[0,m,0]
                        ,[0,0,z]]
    z = konst 0 (p, p)

test_corr2 :: Int -> L.Matrix Float -> L.Matrix Float -> IO (L.Matrix Float)
test_corr2 p k m | w > s     = test_corr2 p m k
                 | otherwise = do x@(U.DenseMatrixArray _ _ _ vx) <- U.newDenseMatrixArray 1 r c
                                  k' <- U.DenseMatrix w h <$> V.thaw (flatten k)
                                  m' <- U.DenseMatrix s t <$> V.thaw (flatten m)
                                  U.corr2 p (BV.singleton k') m' (x U.<<=)
                                  reshape c <$> V.freeze vx
  where
    (w,h) = L.size k
    (s,t) = L.size m
    (r,c) = (s-w+2*p+1, t-h+2*p+1)

test_corr2_arr :: Int -> [L.Matrix Float] -> L.Matrix Float -> IO [L.Matrix Float]
test_corr2_arr p ks m = do x@(U.DenseMatrixArray _ _ _ vx) <- U.newDenseMatrixArray n r c
                           print ("test", n, r, c, MV.length vx)
                           ks' <- mapM (\k -> U.DenseMatrix w h <$> V.thaw (flatten k)) ks
                           m'  <- U.DenseMatrix s t <$> V.thaw (flatten m)
                           U.corr2 p (BV.fromList ks') m' (x U.<<=)
                           let vm = U.denseMatrixArrayToVector x
                           vhm <- BV.mapM (\(U.DenseMatrix _ _ vx) -> reshape c <$> V.freeze vx) vm
                           return $ BV.toList vhm

  where
    n     = length ks
    (w,h) = L.size (head ks)
    (s,t) = L.size m
    (r,c) = (s-w+2*p+1, t-h+2*p+1)

eq :: L.Matrix Float -> L.Matrix Float -> Bool
eq a b = V.all id $ ratio a b

ratio a b =
  let va = flatten a
      vb = flatten b
      ae :: V.Vector Float
      ae = V.zipWith (\a b -> abs (a - b)) va vb
      aa = V.sum ae / fromIntegral (V.length ae)
      err x 0 = x < 0.1
      err x y = let e = x-y
                in (abs (e / y) < 0.02)
  in V.zipWith err va vb

max_stars_repo_name: pierric/neural-network

max_issues_repo_name: pierric/neural-network

max_forks_repo_name: pierric/neural-network


text:
{-# LANGUAGE ViewPatterns #-}
module AsteriskGaussian where

import           Control.Parallel.Strategies
import           Data.Array.Repa                        as R
import           Data.Complex                           as C
import           Data.List                              as L
import           Data.Vector.Storable                   as VS
import           Data.Vector.Unboxed                    as VU
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Math.Gamma
import           Pinwheel.FourierSeries2D
-- import           Pinwheel.Gaussian
import qualified FourierPinwheel.Filtering              as FP
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.BLAS
import           Utils.List
import Utils.Distribution

main = do
  args@(numR2FreqStr:numThetaFreqStr:numRFreqStr:numPointStr:numThetaStr:numRStr:xStr:yStr:periodStr:std1Str:std2Str:stdR2Str:_) <-
    getArgs
  let numR2Freq = read numR2FreqStr :: Int
      numThetaFreq = read numThetaFreqStr :: Int
      numRFreq = read numRFreqStr :: Int
      numPoint = read numPointStr :: Int
      numTheta = read numThetaStr :: Int
      numR = read numRStr :: Int
      x = read xStr :: Double
      y = read yStr :: Double
      period = read periodStr :: Double
      std1 = read std1Str :: Double
      std2 = read std2Str :: Double
      stdR2 = read stdR2Str :: Double
      folderPath = "output/test/AsteriskGaussian"
  let centerR2Freq = div numR2Freq 2
      centerThetaFreq = div numThetaFreq 2
      periodEnv = period ^ 2 / 4
      a = std1
      b = std2
      asteriskGaussianFreqTheta =
        VS.concat .
        parMap
          rdeepseq
          (\tFreq' ->
             let tFreq = tFreq' - centerThetaFreq
                 arr =
                   R.map (* (gaussian1DFreq (fromIntegral tFreq) b :+ 0)) $
                   analyticalFourierCoefficients1
                     numR2Freq
                     1
                     tFreq
                     0
                     a
                     period
                     periodEnv
              in VS.convert . toUnboxed . computeS $
                 centerHollowArray numR2Freq arr) $
        [0 .. numThetaFreq - 1]
      deltaTheta = 2 * pi / Prelude.fromIntegral numTheta
      harmonicsThetaD =
        fromFunction (Z :. numTheta :. numThetaFreq :. numR2Freq :. numR2Freq) $ \(Z :. theta' :. tFreq' :. xFreq' :. yFreq') ->
          let tFreq = fromIntegral $ tFreq' - centerThetaFreq
              xFreq = fromIntegral $ xFreq' - centerR2Freq
              yFreq = fromIntegral $ yFreq' - centerR2Freq
              theta = fromIntegral theta' * deltaTheta
           in (1 / (2 * pi) :+ 0) *
              cis (2 * pi / period * (x * xFreq + y * yFreq) + tFreq * theta)
  harmonicsTheta <- computeUnboxedP harmonicsThetaD
  xs <-
    VS.toList . VS.map magnitude <$>
    gemmBLAS
      numTheta
      1
      (numThetaFreq * numR2Freq ^ 2)
      (VS.convert . toUnboxed $ harmonicsTheta)
      asteriskGaussianFreqTheta
  toFile def (folderPath </> "theta.png") $ do
    layout_title .=
      printf "(%d , %d)" (Prelude.round x :: Int) (Prelude.round y :: Int)
    plot
      (line
         ""
         [ L.zip
             [fromIntegral i * deltaTheta / pi | i <- [0 .. numTheta - 1]]
             xs
         ])
  let centerRFreq = div numRFreq 2
      centerR2Freq = div numR2Freq 2
      gaussian2D =
        computeUnboxedS . fromFunction (Z :. numR2Freq :. numR2Freq) $ \(Z :. i' :. j') ->
          let i = i' - centerR2Freq
              j = j' - centerR2Freq
           in exp
                (pi * fromIntegral (i ^ 2 + j ^ 2) /
                 ((-1) * period ^ 2 * stdR2 ^ 2)) /
              (2 * pi * stdR2 ^ 2) :+
              0
      asteriskGaussianFreqR =
        VS.concat .
        parMap
          rdeepseq
          (\rFreq' ->
             let rFreq = rFreq' - centerRFreq
                 arr =
                   R.map
                     (* ((gaussian1DFourierCoefficients
                            (fromIntegral rFreq)
                            (log periodEnv)
                            b :+
                          0) -- *
                         -- cis
                         --   ((-2) * pi / log periodEnv * fromIntegral rFreq *
                         --    log 0.25)
                        )) . centerHollowArray numR2Freq $
                   analyticalFourierCoefficients1
                     numR2Freq
                     1
                     0
                     rFreq
                     a
                     period
                     periodEnv
              in VS.convert . toUnboxed . computeS . centerHollowArray numR2Freq $
                 arr -- *^ gaussian2D
           ) $
        [0 .. numRFreq - 1]
      deltaR = log periodEnv / Prelude.fromIntegral numR
      harmonicsRD =
        fromFunction (Z :. numR :. numRFreq :. numR2Freq :. numR2Freq) $ \(Z :. r' :. rFreq' :. xFreq' :. yFreq') ->
          let rFreq = fromIntegral $ rFreq' - centerRFreq
              xFreq = fromIntegral $ xFreq' - centerR2Freq
              yFreq = fromIntegral $ yFreq' - centerR2Freq
              r = fromIntegral r' * deltaR - log periodEnv / 2
           in (1 / log periodEnv :+ 0) *
              cis
                (2 * pi / period * (x * xFreq + y * yFreq) +
                 2 * pi / log periodEnv * rFreq * r)
  harmonicsR <- computeUnboxedP harmonicsRD
  ys <-
    VS.toList . VS.map magnitude <$>
    gemmBLAS
      numR
      1
      (numRFreq * numR2Freq ^ 2)
      (VS.convert . toUnboxed $ harmonicsR)
      asteriskGaussianFreqR
  toFile def (folderPath </> "R.png") $ do
    layout_title .= printf "(%.3f , %.3f)" x y
    plot
      (line
         ""
         [ L.zip
             [ fromIntegral i * deltaR - log periodEnv / 2
             | i <- [0 .. numR - 1]
             ]
             ys
         ])


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
-- |
-- Module      : Bench.HMatrix
-- Copyright   : [2017..2020] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Bench.HMatrix ( bench_hmatrix )
  where

import Bench.Util

import Control.DeepSeq
import Criterion.Main
import Data.Proxy
import Foreign.Storable
import Numeric.LinearAlgebra                                        hiding ( randomVector )
import Prelude                                                      hiding ( (<>) )
import System.Random.MWC
import Text.Printf

bench_hmatrix :: [Benchmark]
bench_hmatrix =
  [ bench_level2
  , bench_level3
  ]

bench_level2 :: Benchmark
bench_level2 =
  bgroup "matrix-vector"
    [ bgroup "(#>)"
      [ gemv  200  400
      , gemv  500 1000
      , gemv 1000 2000
      , gemv 2000 3000
      ]
    , bgroup "(<#)"
      [ gevm  200  400
      , gevm  500 1000
      , gevm 1000 2000
      , gevm 2000 3000
      ]
    ]
  where
    gemv :: Int -> Int -> Benchmark
    gemv m n =
      let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e)
          setup _ = withSystemRandom $ \gen -> do
            matA <- randomMatrix gen m n
            vecx <- randomVector gen n
            return (matA, vecx)

          go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark
          go t = env (setup t)
               $ \ ~(matA, vecx) -> bench (showType t)
               $ whnf (matA #>) vecx
      in
      bgroup (printf "%dx%d" m n) (sdcz go)

    gevm :: Int -> Int -> Benchmark
    gevm m n =
      let setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Vector e)
          setup _ = withSystemRandom $ \gen -> do
            matA <- randomMatrix gen m n
            vecx <- randomVector gen m
            return (matA, vecx)

          go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark
          go t = env (setup t)
               $ \ ~(matA, vecx) -> bench (showType t)
               $ whnf (vecx <#) matA
      in
      bgroup (printf "%dx%d" m n) (sdcz go)

bench_level3 :: Benchmark
bench_level3 =
  bgroup "matrix-matrix"
    [ bgroup "(<>)"
      [ gemm  100  100  100
      , gemm  250  250  250
      , gemm  500  500  500
      , gemm 1000 1000 1000
      ]
    ]
  where
    gemm :: Int -> Int -> Int -> Benchmark
    gemm m n k =
      let
          setup :: (Variate e, Storable e) => proxy e -> IO (Matrix e, Matrix e)
          setup _ = withSystemRandom $ \gen -> do
            matA <- randomMatrix gen m k
            matB <- randomMatrix gen k n
            return (matA, matB)

          go :: (Variate e, Numeric e, NFData e, Show (ArgType e)) => proxy e -> Benchmark
          go t = env (setup t)
               $ \ ~(matA, matB) -> bench (showType t)
               $ whnf (matA <>) matB
      in
      bgroup (printf "%dx%dx%d" m n k) (sdcz go)

randomVector :: (Variate e, Storable e) => GenIO -> Int -> IO (Vector e)
randomVector = uniformVector

randomMatrix :: (Variate e, Storable e) => GenIO -> Int -> Int -> IO (Matrix e)
randomMatrix gen m n = do
  v <- uniformVector gen (m * n)
  return $ reshape n v

sdcz :: (forall (e :: *). (Variate e, Numeric e, NFData e, Show (ArgType e)) => Proxy e -> Benchmark)
     -> [Benchmark]
sdcz go =
  [ go (Proxy :: Proxy Float)
  , go (Proxy :: Proxy Double)
  , go (Proxy :: Proxy (Complex Float))
  , go (Proxy :: Proxy (Complex Double))
  ]

max_stars_repo_name: statusfailed/accelerate-blas

max_issues_repo_name: statusfailed/accelerate-blas

max_forks_repo_name: statusfailed/accelerate-blas


text:
{-# OPTIONS_GHC  -fno-warn-unused-binds -fno-warn-unused-matches -fno-warn-name-shadowing -fno-warn-missing-signatures #-}
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses, UndecidableInstances, FlexibleContexts, TypeSynonymInstances #-}

---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------
-- |
-- | Module : Test Morris
-- | Creator: Xiao Ling
-- | Created: 12/8/2015
-- | see    : https://github.com/snoyberg/conduit
-- |
---------------------------------------------------------------------------------------------------
---------------------------------------------------------------------------------------------------

module Tmorris where

import Data.List
import Data.Random
import Data.Conduit
import Data.List.Split
import qualified Data.Conduit.List as Cl

import Control.Monad.Identity

import Criterion.Main

import Core
import Statistics
import Morris

main :: IO ()
main = return ()