text:
-- | https://en.wikibooks.org/wiki/Write_Yourself_a_Scheme_in_48_Hours/Parsing

module Main where

import           Control.Monad
import           Data.Char                     (digitToInt, isDigit)
import           Data.Complex
import           Data.Ratio                    (denominator, numerator, (%))
import           Numeric                       (readDec, readFloat, readHex,
                                                readInt, readOct)
import           System.Environment
import           Text.ParserCombinators.Parsec hiding (spaces)

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

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

-- Whitespace
spaces :: Parser ()
spaces = skipMany1 space

-- Return Values
data LispVal = Atom String
             | List [LispVal]
             | DottedList [LispVal] LispVal
             | Number Integer
             | Float Double
             | Ratio Rational
             | Complex (Complex Double)
             | Character Char
             | String String
             | Bool Bool

parseString :: Parser LispVal
parseString = do
  char '"'
  str <- many (parseHelper <|> noneOf "\\\"")
  char '"'
  return $ String str
  <?> "a string for parseString"
  where parseHelper :: Parser Char
        parseHelper = do
          char '\\'
          x <- oneOf "\\\"nrt"
          return $ case x of
            '\\' -> x
            '"'  -> x
            'n'  -> '\n'
            'r'  -> '\r'
            't'  -> '\t'

parseAtom :: Parser LispVal
parseAtom = do
  first <- letter <|> symbol
  rest <- many $ letter <|> symbol <|> digit
  return $ Atom (first : rest)
  <?> "a atom for parseAtom"

parseBool :: Parser LispVal
parseBool = do
  x <- try (string "#t") <|> string "#f"
  return $ if x == "#t"
           then Bool True
           else Bool False

parseNumber :: Parser LispVal
parseNumber = parseDigit <|> try parseOct <|> try parseDec <|> try parseHex <|> parseBin
  where parseOct = do
          char '#' >> oneOf "oO"
          str <- many1 octDigit
          let [(n, _)] = readOct str
          return $ Number n
        parseDec = do
          char '#' >> oneOf "dD"
          str <- many1 digit
          let [(n, _)] = readDec str
          return $ Number n
        parseHex = do
          char '#' >> oneOf "xX"
          str <- many1 hexDigit
          let [(n, _)] = readHex str
          return $ Number n
        parseBin = do
          char '#' >> oneOf "bB"
          str <- many1 $ oneOf "01"
          let [(n, _)] = readInt 2 (`elem` "01") digitToInt str
          return $ Number n
        parseDigit = do
          str <- many1 digit
          return $ (Number . read) str

parseChar :: Parser LispVal
parseChar = do
  string "#\\"
  x <- many $ letter <|> digit
  return $ case length x of
    0 -> Character '\n'
    1 -> Character $ head x
    _ -> case x of
      "newline" -> Character '\n'
      "tab"     -> Character '\t'
      "space"   -> Character ' '

parseFloat :: Parser LispVal
parseFloat = do
  x <- many1 digit
  char '.'
  y <- many1 digit
  return $ let [(n, _)] = readFloat (x ++ '.' : y) in Float n

parseRatio :: Parser LispVal
parseRatio = do
  x <- many1 digit
  char '/'
  y <- many1 digit
  return $ Ratio (read x % read y)

parseComplex :: Parser LispVal
parseComplex = do
  x <- try parseFloat <|> parseNumber
  char '+'
  y <- try parseFloat <|> parseNumber
  char 'i'
  return $ Complex (toDouble x :+ toDouble y)
  where toDouble :: LispVal -> Double
        toDouble (Float f) = f
        toDouble (Number n) = fromIntegral n

-- Recursive Parsers: Adding lists, dotted lists, and quoted datums

parseList :: Parser LispVal
parseList = do
  char '('
  x <- sepBy parseExpr spaces
  char ')'
  return $ List x

parseDottedList :: Parser LispVal
parseDottedList = do
  char '('
  x <- endBy parseExpr spaces
  y <- char '.' >> spaces >> parseExpr
  char ')'
  return $ DottedList x y

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

parseBackquote :: Parser LispVal
parseBackquote = do
  char '`'
  x <- parseExpr
  return $ List [Atom "quasiquote", x]

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

parseUnquoteSplicing :: Parser LispVal
parseUnquoteSplicing = do
  char ',' >> char '@'
  y <- parseExpr
  return $ List [Atom "unquote-splicing", y]

parseExpr :: Parser LispVal
parseExpr = parseAtom
            <|> try parseBool
            <|> try parseRatio
            <|> try parseComplex
            <|> try parseFloat
            <|> try parseNumber
            <|> parseChar
            <|> parseString
            <|> try parseList
            <|> parseDottedList
            <|> parseQuoted
            <|> parseBackquote
            <|> try parseUnquoteSplicing
            <|> parseUnquote

-- Beginning the Evaluator
showVal :: LispVal -> String
showVal (Atom x) = x
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Ratio x) = show (numerator x) ++ ('/' : show (denominator x))
showVal (Complex x) = show (realPart x) ++ ('+' : show (imagPart x) ++ "i")
showVal (Float x) = show x
showVal (Number x) = show x
showVal (Character x) = "#\\" ++ show x
showVal (String x) = "\"" ++ x ++ "\""
showVal (List x) = "(" ++ unwordsList x ++ ")"
showVal (DottedList first rest)= "(" ++ unwordsList first ++ " . " ++ showVal rest ++ ")"

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

instance Show LispVal where
  show = showVal

readExpr :: String -> LispVal
readExpr input = case parse parseExpr "lisp" input of
  Left err -> String $ "No match: " ++ show err
  Right val -> val

-- Beginnings of an evaluator: Primitives
main :: IO ()
main = getArgs >>= print . eval . readExpr . head

eval :: LispVal -> LispVal
eval val@(String _) = val
eval val@(Bool _) = val
eval val@(Ratio _) = val
eval val@(Complex _) = val
eval val@(Float _) = val
eval val@(Number _) = val
eval val@(Character _) = val
eval (List [Atom "quote", val]) = val

-- Adding basic primitives
eval (List (Atom func : args)) = apply func $ map eval args

apply :: String -> [LispVal] -> LispVal
apply func args = maybe (Bool False) ($ args) $ lookup func primitives

primitives :: [(String, [LispVal] -> LispVal)]
primitives = [("+", numericBinop (+)),
              ("-", numericBinop (-)),
              ("*", numericBinop (*)),
              ("/", numericBinop div),
              ("mod", numericBinop mod),
              ("quotient", numericBinop quot),
              ("remainder", numericBinop rem),
              ("symbol?", unaryOp symbolp),
              ("string?", unaryOp stringp),
              ("number?", unaryOp numberp),
              ("bool?", unaryOp boolp),
              ("list?", unaryOp listp),
              ("pair?", unaryOp pairp)]

numericBinop :: (Integer -> Integer -> Integer) -> [LispVal] -> LispVal
numericBinop op args = Number $ foldl1 op $ map unpackNum args
  where unpackNum :: LispVal -> Integer
        unpackNum (Number n) = n
        unpackNum (String n) = let parsed = readDec n in
          if null parsed
          then 0
          else fst $ head parsed
        unpackNum (List [n]) = unpackNum n
        unpackNum _ = 0

unaryOp :: (LispVal -> LispVal) -> [LispVal] -> LispVal
unaryOp op [x] = op x

symbolp :: LispVal -> LispVal
symbolp (Atom _) = Bool True
symbolp _ = Bool False

stringp :: LispVal -> LispVal
stringp (String _) = Bool True
stringp _ = Bool False

numberp :: LispVal -> LispVal
numberp (Number _) = Bool True
numberp _ = Bool False

boolp :: LispVal -> LispVal
boolp (Bool _) = Bool True
boolp _ = Bool False

listp :: LispVal -> LispVal
listp (List _) = Bool True
listp _ = Bool False

pairp :: LispVal -> LispVal
pairp (List _) = Bool True
pairp (DottedList _ _) = Bool True
pairp _ = Bool False

max_stars_repo_name: zeqing-guo/Haskell-Exercises

max_issues_repo_name: zeqing-guo/Haskell-Exercises

max_forks_repo_name: zeqing-guo/Haskell-Exercises


text:
module Numeric.FFT.Utils
       ( omega, slicevecs, slicemvecs, primes, isPrime
       , allFactors, factors
       , primitiveRoot, invModN, log2, isPow2, dupperm, (%.%)
       , compositions, makeComp, multisetPerms
       , backpermuteM
       ) where

import qualified Control.Monad               as CM
import           Control.Monad.ST
import           Data.Bits
import           Data.Complex
import           Data.List                   (nub)
import qualified Data.List                   as L
import qualified Data.Vector                 as V
import           Data.Vector.Unboxed
import qualified Data.Vector.Unboxed.Mutable as MV
import           Prelude                     hiding (all, concatMap, dropWhile,
                                              enumFromTo, filter, head, length,
                                              map, maximum, null, reverse)
import qualified Prelude                     as P

import           Numeric.FFT.Types

-- | Roots of unity.
omega :: Int -> Complex Double
omega n = cis (2 * pi / fromIntegral n)

-- | Slice a vector @v@ into equally sized parts, each of length @m@.
slicevecs :: Int -> VCD -> VVCD
slicevecs m v = V.map (\i -> slice (i * m) m v) $
                V.enumFromN 0 (length v `div` m)

-- | Slice a mutable vector @v@ into equally sized parts, each of
-- length @m@.
slicemvecs :: Int -> MVCD a -> VMVCD a
slicemvecs m v = V.map (\i -> MV.slice (i * m) m v) $
                 V.enumFromN 0 (MV.length v `div` m)

-- | Determine primitive roots modulo n.
--
-- From Wikipedia (https://en.wikipedia.org/wiki/Primitive_root_modulo_n):
--
-- No simple general formula to compute primitive roots modulo n is
-- known.  There are however methods to locate a primitive root that
-- are faster than simply trying out all candidates.  If the
-- multiplicative order of a number m modulo n is equal to phi(n) (the
-- order of Z_n^x), then it is a primitive root.  In fact the converse
-- is true: if m is a primitive root modulo n, then the multiplicative
-- order of m is phi(n).  We can use this to test for primitive roots.
-- [Here, phi(n) is Euler's totient function, and Z_n^x is the
-- multiplicative group of integers modulo n.]
--
-- First, compute phi(n).  Then determine the different prime factors
-- of phi(n), say p1, ..., pk.  Now, for every element m of Z_n^x,
-- compute
--
--    m^(phi(n) / pi) mod n   for i = 1, ..., k
--
-- using a fast algorithm for modular exponentiation such as
-- exponentiation by squaring.  A number m for which these k results
-- are all different from 1 is a primitive root.
--
-- [In our case, n is restricted to being prime, and phi(p) = p - 1
-- for prime p.]
--
primitiveRoot :: Int -> Int
primitiveRoot p
  | isPrime p =
    let tot = p - 1
        -- ^ Euler's totient function for prime values.
        totpows = map (tot `div`) $ fromList $ nub $ toList $ allFactors tot
        -- ^ Powers to check.
        check n = all (/=1) $ map (expt p n) totpows
        -- ^ All powers are different from 1 => primitive root.
    in fromIntegral $ head $ dropWhile (not . check) $ fromList [1..p-1]
  | otherwise = error "Attempt to take primitive root of non-prime value"

-- | Fast exponentation modulo n by squaring.
expt :: Int -> Int -> Int -> Int
expt n b pow = fromIntegral $ go pow
  where bb = fromIntegral b
        nb = fromIntegral n
        go :: Int -> Integer
        go p
          | p == 0 = 1
          | p `mod` 2 == 1 = (bb * go (p - 1)) `mod` nb
          | otherwise = let h = go (p `div` 2) in (h * h) `mod` nb

-- | Find inverse element in multiplicative integer group modulo n.
invModN :: Int -> Int -> Int
invModN n g = head $ filter (\iv -> (g * iv) `mod` n == 1) $ enumFromTo 1 (n-1)

-- | Prime sieve from Haskell wiki.
primes :: Integral a => [a]
primes = 2 : primes'
  where primes' = sieve [3, 5 ..] 9 primes'
        sieve (x:xs) q ps@ ~(p:t)
          | x < q     = x : sieve xs q ps
          | otherwise =     sieve [n | n <- xs, rem n p /= 0] (P.head t^(2::Int)) t

-- | Naive primality testing.
isPrime :: Integral a => a -> Bool
isPrime n = n `P.elem` P.takeWhile (<= n) primes

-- | Simple prime factorisation.
allFactors :: (Integral a, Unbox a) => a -> Vector a
allFactors n = fromList $ go n primes
  where go cur pss@(p:ps)
          | cur == p         = [p]
          | cur `mod` p == 0 = p : go (cur `div` p) pss
          | otherwise        = go cur ps

-- | Simple prime factorisation: small factors only; largest/last
-- factor picked out as "special".
factors :: (Integral a, Unbox a) => a -> (a, Vector a)
factors n = let (lst, rest) = go n primes in (lst, fromList rest)
  where go cur pss@(p:ps)
          | cur == p         = (p, [])
          | cur `mod` p == 0 = let (lst, rest) = go (cur `div` p) pss
                               in (lst, p : rest)
          | otherwise        = go cur ps

-- | Base-2 logarithm.
log2 :: Int -> Int
log2 1 = 0
log2 n = 1 + log2 (n `div` 2)

-- | Check for powers of two.
isPow2 :: Int -> Bool
isPow2 1 = True
isPow2 n
  | n `mod` 2 == 0 = isPow2 $ n `div` 2
  | otherwise      = False

-- | Duplicate a sub-permutation to fill a given vector length.
dupperm :: Int -> VI -> VI
dupperm n p =
  let sublen = length p
      shift di = map (+(sublen * di)) p
  in concatMap shift $ enumFromN 0 (n `div` sublen)

-- | Composition of permutations.
(%.%) :: VI -> VI -> VI
p1 %.% p2 = backpermute p2 p1

-- | Generate all compositions of a given integer.
compositions :: Int -> V.Vector (Vector Int)
compositions 0 = V.empty
compositions n = let fs = allFactors n
                 in V.reverse $ V.map (makeComp fs) $ V.enumFromN 0 (2^(n-1))

-- | Generate a single composition of a given integer.
makeComp :: Vector Int -> Int -> Vector Int
makeComp fs i = fromList $ foldOps (toList fs) $ makeOps (length fs) i
  where foldOps :: [Int] -> [Bool] -> [Int]
        foldOps (f:fs) ops = go f fs ops
          where go acc [] [] = [acc]
                go acc (f:fs) (op:ops) = if op then go (acc * f) fs ops
                                         else acc : go f fs ops
        makeOps :: Int -> Int -> [Bool]
        makeOps n i = P.replicate (n - 1 - P.length bs) False P.++ bs
          where bs = P.dropWhile not $ P.reverse $
                     P.map (testBit i) [0..finiteBitSize i-1]

-- | Generate all distinct permutations of a multiset in lexicographic
-- order.
multisetPerms :: Vector Int -> [Vector Int]
multisetPerms idp = sidp : L.unfoldr step sidp
  where sidp = fromList $ L.sort $ toList idp
        step v = case permStep v of
          Nothing -> Nothing
          Just p  -> Just (p, p)
        permStep :: Vector Int -> Maybe (Vector Int)
        permStep v =
          if null ks
          then Nothing
          else let k = maximum ks
                   ls = filter (\i -> v ! k < v ! i) $ enumFromN 0 n
                   l = maximum ls
               in Just $ revEnd k (swap k l)
          where n = length v
                ks = filter (\i -> v ! i < v ! (i+1)) $ enumFromN 0 (n-1)
                swap a b = generate n $ \i ->
                  if i == a then v ! b
                  else if i == b then v ! a
                       else v ! i
                revEnd f vv = generate n $ \i ->
                  if i <= f
                  then vv ! i
                  else vv ! (n - i + f)

-- | ST monad version of vector permutation.
backpermuteM :: Int -> VI -> MVCD s -> MVCD s -> ST s ()
backpermuteM n perm vin vout =
  CM.forM_ [0..n-1] $ \i -> do
    idx <- indexM perm i
    x <- MV.unsafeRead vin idx
    MV.unsafeWrite vout i x

ext: hs

max_stars_repo_name: ian-ross/arb-fft

max_issues_repo_name: ian-ross/arb-fft

max_forks_repo_name: ian-ross/arb-fft


text:
{-# language InstanceSigs, DerivingStrategies #-}
{-# language PartialTypeSignatures #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures #-}
{-| This module defines the 'Shaped' typeclass, which is used to generically
    manipulate values as fixed-points of higher-order functors in order to
    analyze their structure, e.g. while observing evaluation.

    If you just care about testing the strictness of functions over datatypes
    which are already instances of @Shaped@, you don't need to use this module.

    __Important note:__ To define new instances of 'Shaped' for types which
    implement 'GHC.Generic', __an empty instance will suffice__, as all the
    methods of 'Shaped' can be filled in by generic implementations. For
    example:

    > import GHC.Generics as GHC
    > import Generics.SOP as SOP
    >
    > data D = C deriving (GHC.Generic)
    >
    > instance SOP.Generic D
    > instance SOP.HasDatatypeInfo D
    >
    > instance Shaped D

    Using the @DeriveAnyClass@ extension, this can be shortened to one line:

    > data D = C deriving (GHC.Generic, SOP.Generic, SOP.HasDatatypeInfo, Shaped)

    Manual instances of 'Shaped' are necessary for types which do not or cannot
    implement GHC's @Generic@ typeclass, such as existential types, abstract
    types, and GADTs.

    This module is heavily based upon the approach in "Data.Functor.Foldable",
    which in turn is modeled after the paper "Functional Programming with
    Bananas, Lenses, Envelopes and Barbed Wire" (1991) by Erik Meijer, Maarten
    Fokkinga and Ross Paterson. If you don't yet understand recursion schemes
    and want to understand this module, it's probably a good idea to familiarize
    yourself with "Data.Functor.Foldable" before diving into this higher-order
    generalization.
-}
module Test.StrictCheck.Shaped
  ( Shaped(..)
  , module Test.StrictCheck.Shaped.Flattened
  -- * Fixed-points of 'Shape's
  , type (%)(..)
  -- * Folds and unfolds over fixed-points of @Shape@s
  , unwrap
  , interleave
  , (%)
  , fuse
  , translate
  , fold
  , unfold
  , unzipWith
  -- , reshape
  -- * Rendering 'Shaped' things as structured text
  , QName
  , Rendered(..)
  , RenderLevel(..)
  , renderfold
  -- * Tools for manually writing instances of 'Shaped'
  -- ** Implementing 'Shaped' for primitive types
  , Prim(..), unPrim
  , projectPrim
  , embedPrim
  , matchPrim
  , flatPrim
  , renderPrim
  , renderConstant
  -- ** Implementing 'Shaped' for container types
  , Containing(..)
  , projectContainer
  , embedContainer
  -- * Generic implementation of the methods of 'Shaped'
  , GShaped
  , GShape(..)
  , gProject
  , gEmbed
  , gMatch
  , gRender
  ) where

import Type.Reflection
import Data.Functor.Product
import Data.Bifunctor
import Data.Bifunctor.Flip
import Data.Coerce

import Generics.SOP hiding ( Shape )

import Data.Complex
-- import Data.List.NonEmpty (NonEmpty(..))

import Test.StrictCheck.Shaped.Flattened

-- TODO: provide instances for all of Base

-- | When a type @a@ is @Shaped@, we know how to convert it into a
-- representation parameterized by an arbitrary functor @f@, so that @Shape a f@
-- (the "shape of @a@ parameterized by @f@") is structurally identical to the
-- topmost structure of @a@, but with @f@ wrapped around any subfields of @a@.
--
-- Note that this is /not/ a recursive representation! The functor @f@ in
-- question wraps the original type of the field and /not/ a @Shape@ of that
-- field.
--
-- For instance, the @Shape@ of @Either a b@ might be:
--
-- > data EitherShape a b f
-- >   = LeftShape  (f a)
-- >   | RightShape (f b)
-- >
-- > instance Shaped (Either a b) where
-- >   type Shape (Either a b) = EitherShape a b
-- >   ...
--
-- The shape of a primitive type should be isomorphic to the primitive type,
-- with the functor parameter left unused.
class Typeable a => Shaped (a :: *) where
  -- | The @Shape@ of an @a@ is a type isomorphic to the outermost level of
  -- structure in an @a@, parameterized by the functor @f@, which is wrapped
  -- around any fields (of any type) in the original @a@.
  type Shape a :: (* -> *) -> *
  type Shape a = GShape a

  -- | Given a function to expand any @Shaped@ @x@ into an @f x@, expand an @a@
  -- into a @Shape a f@
  --
  -- That is: convert the top-most level of structure in the given @a@ into a
  -- @Shape@, calling the provided function on each field in the @a@ to produce
  -- the @f x@ necessary to fill that hole in the produced @Shape a f@.
  --
  -- Inverse to 'embed'.
  project :: (forall x. Shaped x => x -> f x) -> a -> Shape a f

  default project
    :: GShaped a
    => (forall x. Shaped x => x -> f x)
    -> a
    -> Shape a f
  project = gProject

  -- | Given a function to collapse any @f x@ into a @Shaped@ @x@, collapse a
  -- @Shape a f@ into merely an @a@
  --
  -- That is: eliminate the top-most @Shape@ by calling the provided function on
  -- each field in that @Shape a f@, and using the results to fill in the pieces
  -- necessary to build an @a@.
  --
  -- Inverse to 'project'.
  embed :: (forall x. Shaped x => f x -> x) -> Shape a f -> a

  default embed
    :: GShaped a
    => (forall x. Shaped x => f x -> x)
    -> Shape a f
    -> a
  embed = gEmbed

  -- | Given two @Shape@s of the same type @a@ but parameterized by potentially
  -- different functors @f@ and @g@, pattern-match on them to expose a uniform
  -- view on their fields (a 'Flattened' @(Shape a)@) to a continuation which
  -- may operate on those fields to produce some result
  --
  -- If the two supplied @Shape@s do not structurally match, only the fields of
  -- the first are given to the continuation. If they do match, the fields of
  -- the second are also given, along with type-level proof that the types of
  -- the two sets of fields align.
  --
  -- This very general operation subsumes equality testing, mapping, zipping,
  -- shrinking, and many other structural operations over @Shaped@ things.
  --
  -- It is somewhat difficult to manually write instances for this method, but
  -- consulting its generic implementation 'gMatch' may prove helpful.
  --
  -- See "Test.StrictCheck.Shaped.Flattened" for more information.
  match :: Shape a f -> Shape a g
        -> (forall xs. All Shaped xs
              => Flattened (Shape a) f xs
              -> Maybe (Flattened (Shape a) g xs)
              -> result)
        -> result

  default match :: GShaped a
        => Shape a f -> Shape a g
        -> (forall xs. All Shaped xs
              => Flattened (Shape a) f xs
              -> Maybe (Flattened (Shape a) g xs)
              -> result)
        -> result
  match = gMatch

  -- | Convert a @Shape a@ whose fields are some unknown constant type into a
  -- 'RenderLevel' filled with that type
  --
  -- This is a specialized pretty-printing mechanism which allows for displaying
  -- counterexamples in a structured format. See the documentation for
  -- 'RenderLevel'.
  render :: Shape a (K x) -> RenderLevel x

  default render :: (GShaped a, HasDatatypeInfo a)
          => Shape a (K x) -> RenderLevel x
  render = gRender

-- * Fixed-points of 'Shape's

-- | A value of type @f % a@ has the same structure as an @a@, but with the
-- structure of the functor @f@ interleaved at every field (including ones of
-- types other than @a@). Read this type aloud as "a interleaved with f's".
newtype (f :: * -> *) % (a :: *) :: * where
  Wrap :: f (Shape a ((%) f)) -> f % a

-- | Look inside a single level of an interleaved @f % a@. Inverse to the 'Wrap'
-- constructor.
unwrap :: f % a -> f (Shape a ((%) f))
unwrap (Wrap fs) = fs

-- * Folds and unfolds over fixed-points of @Shape@s

-- | Map a function across all the fields in a 'Shape'
--
-- This function may change the functor over which the @Shape@ is parameterized.
-- It can assume recursively that all the fields in the @Shape@ are themselves
-- instances of @Shaped@ (which they should be!). This means that you can nest
-- calls to @translate@ recursively.
translate :: forall a f g. Shaped a
          => (forall x. Shaped x => f x -> g x)
          -> Shape a f -> Shape a g
translate t d = match @a d d $ \flat _ ->
  unflatten $ mapFlattened @Shaped t flat

-- | The equivalent of a fold (catamorphism) over recursively 'Shaped' values
--
-- Given a function which folds an @f@ containing some @Shape x g@ into a @g x@,
-- recursively fold any interleaved @f % a@ into a @g a@.
fold :: forall a f g. (Functor f, Shaped a)
     => (forall x. Shaped x => f (Shape x g) -> g x)
     -> f % a -> g a
fold alg = alg . fmap (translate @a (fold alg)) . unwrap

-- | The equivalent of an unfold (anamorphism) over recursively 'Shaped' values
--
-- Given a function which unfolds an @f x@ into a @g@ containing some @Shape x
-- f@, corecursively unfold any @f a@ into an interleaved @g % a@.
unfold :: forall a f g. (Functor g, Shaped a)
       => (forall x. Shaped x => f x -> g (Shape x f))
       -> f a -> g % a
unfold coalg = Wrap . fmap (translate @a (unfold coalg)) . coalg

-- TODO: mapM, foldM, unfoldM, ...

-- | Fuse the interleaved @f@-structure out of a recursively interleaved @f %
-- a@, given some way of fusing a single level @f x -> x@.
--
-- This is a special case of 'fold'.
fuse
  :: (Functor f, Shaped a)
  => (forall x. f x -> x)
  -> (f % a -> a)
fuse e = e . fold (fmap (embed e))

-- | Interleave an @f@-structure at every recursive level of some @a@, given
-- some way of generating a single level of structure @x -> f x@.
--
-- This is a special case of 'unfold'.
interleave
  :: (Functor f, Shaped a)
  => (forall x. x -> f x)
  -> (a -> f % a)
interleave p = unfold (fmap (project p)) . p

-- | An infix synonym for 'interleave'
(%) :: forall a f. (Functor f, Shaped a)
    => (forall x. x -> f x)
    -> a -> f % a
(%) = interleave

-- | A higher-kinded @unzipWith@, operating over interleaved structures
--
-- Given a function splitting some @f x@ into a functor-product @Product g h x@,
-- recursively split an interleaved @f % a@ into two interleaved structures:
-- one built of @g@-shapes and one of @h@-shapes.
--
-- Note that @Product ((%) g) ((%) h) a@ is isomorphic to @(g % a, h % a)@; to
-- get the latter, pattern-match on the 'Pair' constructor of 'Product'.
unzipWith
  :: (All Functor [f, g, h], Shaped a)
  => (forall x. f x -> (g x, h x))
  -> (f % a -> (g % a, h % a))
unzipWith split =
  unPair . fold (crunch . pair . split)
  where
    crunch
      :: forall x g h.
      (Shaped x, Functor g, Functor h)
      => Product g h (Shape x (Product ((%) g) ((%) h)))
      -> Product ((%) g) ((%) h) x
    crunch =
      pair
      . bimap (Wrap . fmap (translate @x (fst . unPair)))
              (Wrap . fmap (translate @x (snd . unPair)))
      . unPair

    pair :: (l x, r x) -> Product l r x
    pair = uncurry Pair

    unPair :: Product l r x -> (l x, r x)
    unPair (Pair lx rx) = (lx, rx)

-- | TODO: document this strange function
{-
reshape :: forall b a f g. (Shaped a, Shaped b, Functor f)
        => (f (Shape b ((%) g)) -> g (Shape b ((%) g)))
        -> (forall x. Shaped x => f % x -> g % x)
        -> f % a -> g % a
reshape homo hetero d =
  case eqTypeRep (typeRep @a) (typeRep @b) of
    Nothing    -> hetero d
    Just HRefl ->
      Wrap
      $ homo . fmap (translate @a (reshape @b homo hetero))
      $ unwrap d
-}

----------------------------------
-- Rendering shapes for display --
----------------------------------

-- | Convert an @f % a@ into a structured pretty-printing representation,
-- suitable for further display/processing
renderfold
  :: forall a f. (Shaped a, Functor f)
  => f % a -> Rendered f
renderfold = unK . fold oneLevel
  where
    oneLevel :: forall x. Shaped x
             => f (Shape x (K (Rendered f)))
             -> K (Rendered f) x
    oneLevel = K . RWrap . fmap (render @x)

-- | A @QName@ is a qualified name
--
-- Note:
-- > type ModuleName   = String
-- > type DatatypeName = String
type QName = (ModuleName, DatatypeName, String)

-- | @RenderLevel@ is a functor whose outer shape contains all the information
-- about how to pretty-format the outermost @Shape@ of some value. We use
-- parametricity to make it difficult to construct incorrect 'render' methods,
-- by asking the user merely to produce a single @RenderLevel@ and stitching
-- nested @RenderLevel@s into complete 'Rendered' trees.
data RenderLevel x
  = ConstructorD QName [x]
  -- ^ A prefix constructor, and a list of its fields
  | InfixD QName Associativity Fixity x x
  -- ^ An infix constructor, its associativity and fixity, and its two fields
  | RecordD QName [(QName, x)]
  -- ^ A record constructor, and a list of its field names paired with fields
  | CustomD Fixity
    [Either (Either String (ModuleName, String)) (Fixity, x)]
  -- ^ A custom pretty-printing representation (i.e. for abstract types), which
  -- records a fixity and a list of tokens of three varieties: 1) raw strings,
  -- 2) qualified strings (from some module), or 3) actual fields, annotated
  -- with their fixity
  deriving (Eq, Ord, Show, Functor)

-- | @Rendered f@ is the fixed-point of @f@ composed with 'RenderLevel': it
-- alternates between @f@ shapes and @RenderLevel@s. Usually, @f@ will be the
-- identity functor 'I', but not always.
data Rendered f
  = RWrap (f (RenderLevel (Rendered f)))

----------------------------------------------------
-- Tools for manually writing instances of Shaped --
----------------------------------------------------

-- | The @Shape@ of a spine-strict container (i.e. a @Map@ or @Set@) is the same
-- as a container of demands on its elements. However, this does not have the
-- right /kind/ to be used as a @Shape@.
--
-- The @Containing@ newtype solves this problem. By defining the @Shape@ of some
-- container @(C a)@ to be @(C `Containing` a)@, you can use the methods
-- @projectContainer@ and @embedContainer@ to implement @project@ and @embed@
-- for your container type (although you will still need to manually define
-- @match@ and @render@).
newtype Containing h a f
  = Container (h (f a))
  deriving (Eq, Ord, Show)

-- | Generic implementation of @project@ for any container type whose @Shape@
-- is represented as a @Containing@ newtype
projectContainer :: (Functor c, Shaped a)
  => (forall x. Shaped x => x -> f x)
  -> c a -> Containing c a f
projectContainer p x = Container (fmap p x)

-- | Generic implementation of @embed@ for any container type whose @Shape@
-- is represented as a @Containing@ newtype
embedContainer :: (Functor c, Shaped a)
  => (forall x. Shaped x => f x -> x)
  -> Containing c a f -> c a
embedContainer e (Container x) = fmap e x

-- TODO: helper functions for matching and prettying containers

-- | The @Shape@ of a primitive type should be equivalent to the type itself.
-- However, this does not have the right /kind/ to be used as a @Shape@.
--
-- The @Prim@ newtype solves this problem. By defining the @Shape@ of some
-- primitive type @p@ to be @Prim p@, you can use the methods @projectPrim@,
-- @embedPrim@, @matchPrim@, and @prettyPrim@ to completely fill in the
-- definition of the @Shaped@ class for a primitive type.
--
-- __Note:__ It is only appropriate to use this @Shape@ representation when a
-- type really is primitive, in that it contains no interesting substructure.
-- If you use the @Prim@ representation inappropriately, StrictCheck will not be
-- able to inspect the richer structure of the type in question.
newtype Prim (x :: *) (f :: * -> *)
  = Prim x
  deriving (Eq, Ord, Show)
  deriving newtype (Num)

-- | Get the wrapped @x@ out of a @Prim x f@ (inverse to the @Prim@ constructor)
unPrim :: Prim x f -> x
unPrim (Prim x) = x

-- | Generic implementation of @project@ for any primitive type whose @Shape@ is
-- is represented as a @Prim@ newtype
projectPrim :: (forall x. Shaped x => x -> f x) -> a -> Prim a f
projectPrim _ = Prim

-- | Generic implementation of @embed@ for any primitive type whose @Shape@ is
-- is represented as a @Prim@ newtype
embedPrim :: (forall x. Shaped x => f x -> x) -> Prim a f -> a
embedPrim _ = unPrim

-- | Generic implementation of @match@ for any primitive type whose @Shape@ is
-- is represented as a @Prim@ newtype with an underlying @Eq@ instance
matchPrim :: Eq a => Prim a f -> Prim a g
           -> (forall xs. All Shaped xs
                => Flattened (Prim a) f xs
                -> Maybe (Flattened (Prim a) g xs)
                -> result)
           -> result
matchPrim (Prim a) (Prim b) k =
  k (flatPrim a)
     (if a == b then (Just (flatPrim b)) else Nothing)

-- | Helper for writing @match@ instances for primitive types which don't have
-- @Eq@ instance
--
-- This generates a @Flattened@ appropriate for using in the implementation of
-- @match@. For more documentation on how to use this, see the documentation of
-- 'match'.
flatPrim :: a -> Flattened (Prim a) g '[]
flatPrim x = Flattened (const (Prim x)) Nil

-- | Generic implementation of @render@ for any primitive type whose @Shape@ is
-- is represented as a @Prim@ newtype
renderPrim :: Show a => Prim a (K x) -> RenderLevel x
renderPrim (Prim a) = renderConstant (show a)

-- | Given some @string@, generate a custom pretty-printing representation which
-- just shows the string
renderConstant :: String -> RenderLevel x
renderConstant s = CustomD 11 [Left (Left s)]

-- TODO: What about demands for abstract types with > 1 type of unbounded-count field?

{-
withFieldsContainer ::
  forall c a f result.
     (forall r h.
        c (h a) ->
        (forall x. Shaped x
           => [h x]
           -> (forall g. [g x] -> c (g a))
           -> r)
        -> r)
  -> Containing c a f
  -> (forall xs. All Shaped xs
        => NP f xs
        -> (forall g. NP g xs -> Containing c a g)
        -> result)
  -> result
withFieldsContainer viaContaining (Container c) cont =
  viaContaining c $
    \list un ->
       withNP @Shaped list (Container . un) cont

-- TODO: Make this work for any number of lists of fields, by carefully using
-- unsafeCoerce to deal with unknown list lengths

withFieldsViaList ::
  forall demand f result.
     (forall r h.
        demand h ->
        (forall x. Shaped x
           => [h x]
           -> (forall g. [g x] -> demand g)
           -> r)
        -> r)
  -> demand f
  -> (forall xs. All Shaped xs
        => NP f xs
        -> (forall g. NP g xs -> demand g)
        -> result)
  -> result
withFieldsViaList viaList demand cont =
  viaList demand $
    \list un ->
       withNP @Shaped list un cont

withNP :: forall c demand result f x. c x
       => [f x]
       -> (forall g. [g x] -> demand g)
       -> (forall xs. All c xs
             => NP f xs -> (forall g. NP g xs -> demand g) -> result)
       -> result
withNP list unList cont =
  withUnhomogenized @c list $ \np ->
    cont np (unList . homogenize)

withConcatenated :: NP (NP f) xss -> (forall xs. NP f xs -> r) -> r
withConcatenated pop cont =
  case pop of
    Nil         -> cont Nil
    (xs :* xss) -> withConcatenated xss (withPrepended xs cont)
  where
    withPrepended ::
      NP f ys -> (forall zs. NP f zs -> r)
              -> (forall zs. NP f zs -> r)
    withPrepended pre k rest =
      case pre of
        Nil        -> k rest
        (x :* xs)  -> withPrepended xs (k . (x :*)) rest

homogenize :: All ((~) a) as => NP f as -> [f a]
homogenize      Nil  = []
homogenize (a :* as) = a : homogenize as

withUnhomogenized :: forall c a f r.
  c a => [f a] -> (forall as. (All c as, All ((~) a) as) => NP f as -> r) -> r
withUnhomogenized      []  k = k Nil
withUnhomogenized (a : as) k =
  withUnhomogenized @c as $ \np -> k (a :* np)
-}

---------------------------------------------------
-- Generic implementation of the Shaped methods --
---------------------------------------------------

-- | The 'Shape' used for generic implementations of 'Shaped'
--
-- This wraps a sum-of-products representation from "Generics.SOP".
newtype GShape a f
  = GS (NS (NP f) (Code a))

-- | The collection of constraints necessary for a type to be given a generic
-- implementation of the 'Shaped' methods
type GShaped a =
  ( Generic a
  , Shape a ~ GShape a
  , All2 Shaped (Code a)
  , SListI (Code a)
  , All SListI (Code a) )

-- | Generic 'project'
gProject :: GShaped a
         => (forall x. Shaped x => x -> f x)
         -> a -> Shape a f
gProject p !(from -> sop) =
  GS (unSOP (hcliftA (Proxy @Shaped) (p . unI) sop))

-- | Generic 'embed'
gEmbed :: GShaped a
       => (forall x. Shaped x => f x -> x)
       -> Shape a f -> a
gEmbed e !(GS d) =
  to (hcliftA (Proxy @Shaped) (I . e) (SOP d))

-- | Generic 'match'
gMatch :: forall a f g result. GShaped a
       => Shape a f -> Shape a g
       -> (forall xs. All Shaped xs
             => Flattened (Shape a) f xs
             -> Maybe (Flattened (Shape a) g xs)
             -> result)
       -> result
gMatch !(GS df) !(GS dg) cont =
  go @(Code a) df (Just dg) $ \flatF mflatG ->
    cont (flatGD flatF) (flatGD <$> mflatG)
  where
    go :: forall xss r.
      (All SListI xss, All2 Shaped xss)
       => NS (NP f) xss
       -> Maybe (NS (NP g) xss)
       -> (forall xs. All Shaped xs
             => Flattened (Flip SOP xss) f xs
             -> Maybe (Flattened (Flip SOP xss) g xs)
             -> r)
       -> r
    go (Z (fieldsF :: _ xs)) (Just (Z fieldsG)) k =
      k @xs (flatZ fieldsF)  (Just (flatZ fieldsG))
    go (Z (fieldsF :: _ xs)) _ k =   -- Nothing | Just (S _)
      k @xs (flatZ fieldsF)  Nothing
    go (S moreF) Nothing k =
      go moreF Nothing $ \(flatF :: _ xs) _ ->
        k @xs (flatS flatF) Nothing
    go (S moreF) (Just (Z _)) k =
      go moreF Nothing $ \(flatF :: _ xs) _ ->
        k @xs (flatS flatF) Nothing
    go (S moreF) (Just (S moreG)) k =
      go moreF (Just moreG) $ \(flatF :: _ xs) mflatG ->
        k @xs (flatS flatF) (flatS <$> mflatG)

    flatZ
      :: forall h xs xss. NP h xs -> Flattened (Flip SOP (xs : xss)) h xs
    flatZ = Flattened (Flip . SOP . Z)

    flatS
      :: forall h xs xs' xss.
      Flattened (Flip SOP xss) h xs
      -> Flattened (Flip SOP (xs' : xss)) h xs
    flatS (Flattened un fields) =
      Flattened (Flip . SOP . S . coerce . un) fields

    flatGD :: forall t h xs.
      Flattened (Flip SOP (Code t)) h xs -> Flattened (GShape t) h xs
    flatGD (Flattened un fields) =
      Flattened (GS . coerce . un) fields

-- | Generic 'render'
gRender :: forall a x. (HasDatatypeInfo a, GShaped a)
         => Shape a (K x) -> RenderLevel x
gRender (GS demand) =
  case info of
    ADT m d cs ->
      renderC m d demand cs
    Newtype m d c ->
      renderC m d demand (c :* Nil)
  where
    info = datatypeInfo (Proxy @a)

    renderC :: forall as. ModuleName -> DatatypeName
            -> NS (NP (K x)) as
            -> NP ConstructorInfo as
            -> RenderLevel x
    renderC m d subShape constructors =
      case (subShape, constructors) of
        (Z demandFields, c :* _) ->
          case c of
            Constructor name ->
              ConstructorD (m, d, name) $
                hcollapse demandFields
            Infix name associativity fixity ->
              case demandFields of
                (K a :* K b :* Nil) ->
                  InfixD (m, d, name) associativity fixity a b
            Record name fieldsInfo ->
              RecordD (m, d, name) $
                zip ( hcollapse
                    . hliftA (\(FieldInfo f) -> K (m, d, f))
                    $ fieldsInfo )
                    (hcollapse demandFields)
        (S another, _ :* different) ->
          renderC m d another different

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

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

instance (Typeable a, Typeable b) => Shaped (a -> b) where
  type Shape (a -> b) = Prim (a -> b)
  project = projectPrim
  embed = embedPrim
  match (Prim f) (Prim g) k = k (flatPrim f) (Just $ flatPrim g)
  render _ = renderConstant ("<function> :: " ++ show (typeRep @(a -> b)))

instance Shaped Char where
  type Shape Char = Prim Char
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Word where
  type Shape Word = Prim Word
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Int where
  type Shape Int = Prim Int
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Double where
  type Shape Double = Prim Double
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Float where
  type Shape Float = Prim Float
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Rational where
  type Shape Rational = Prim Rational
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance Shaped Integer where
  type Shape Integer = Prim Integer
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

instance (Typeable a, Eq a, Show a) => Shaped (Complex a) where
  type Shape (Complex a) = Prim (Complex a)
  project = projectPrim
  embed   = embedPrim
  match   = matchPrim
  render  = renderPrim

-- instance Generic (NonEmpty a)
-- instance HasDatatypeInfo (NonEmpty a)
-- instance Shaped a => Shaped (NonEmpty a) where

-- Tree
-- Map k
-- Seq
-- Set
-- IntMap
-- IntSet

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

max_stars_repo_name: janpath/StrictCheck

max_issues_repo_name: janpath/StrictCheck

max_forks_repo_name: janpath/StrictCheck


text:
{-# OPTIONS_GHC -fglasgow-exts #-}
{-# LANGUAGE UndecidableInstances #-}
{- OPTIONS_GHC -fwarn-incomplete-patterns -}

module GeneralizedSignals where

--import
import qualified Data.StorableVector as SV
import qualified Data.StorableVector.Base as SVB
import qualified Data.StorableVector.ST.Strict as SVST
import Foreign.Storable
--import EvalM
import Foreign.Storable.Tuple
import Data.List
import Control.Monad.ST
import Control.Monad
import TNUtils
import VectorsL
import Nats
--import PlotGnuplot
import Data.Ix
import Data.Complex
import Numeric.FFT
import Math.FFT
import qualified Data.Array.CArray as CA
import System.IO.Unsafe

newtype Time = Time Double
    deriving (Eq, Show, Num, Fractional, Floating, Discretizable, IsDouble, Storable)
newtype Freq = Freq Double
    deriving (Eq, Show, Num, Fractional, Floating, Discretizable, IsDouble, Storable)
newtype Length = Length Double
    deriving (Eq, Show, Num, Fractional, Floating, Discretizable, IsDouble, Storable)
newtype SpaceFreq = SpaceFreq Double
    deriving (Eq, Show, Num, Fractional, Floating, Discretizable, IsDouble, Storable)

class (Eq a,Show a,Num a,Fractional a,Floating a)=> IsDouble a where
    toDouble :: a -> Double
    fromDouble :: Double -> a

instance IsDouble Double where
    toDouble = id
    fromDouble = id

class Discretizable a where
    discreteIndex ::  (a,a) -> a -> a -> Int
    discreteRange :: (a,a) -> a -> [a]
    discreteRngLength :: (a,a) -> a -> Int
    foldRangeM :: Monad m => (a,a) -> a -> b -> (a -> b -> m b) -> m b

instance Discretizable Double where
    discreteIndex (lo,hi) dx x = round $ (x-lo)/dx
    discreteRange (lo,hi) dx = [lo, lo+dx..hi]
    discreteRngLength (lo,hi) dx = round $ (hi-lo)/dx
    foldRangeM (lo,hi) dx init f =
        let nmax = discreteRngLength (lo,hi) dx
            go n acc | n > nmax = return acc
                     | otherwise = do
              next <- f ((realToFrac n)*dx+lo) acc
              go (n+1) next
        in go 0 init

instance Discretizable Int where
    discreteIndex (lo,hi) 1 x = x-lo
    discreteIndex (lo,hi) dx x = (x-lo) `div` dx
    discreteRange (lo,hi) dx = [lo,lo+dx..hi]
    discreteRngLength (lo,hi) dx = (hi-lo) `div` dx
    foldRangeM (lo,hi) dx init f =
        let go n acc | n > hi = return acc
                     | otherwise = do
              next <- f n acc
              go (n+dx) next
        in go lo init

{-instance Discretizable a => Discretizable (Vec Z a) where
    discreteIndex (lo,hi) dx x = 0
    discreteRange (lo,hi) dx = [vnil]
    discreteRngLength (lo,hi) dx = 1
    foldRangeM (lo,hi) dx init f = return init -}

instance Discretizable a => Discretizable (Vec (S Z) a) where
    discreteIndex (lo,hi) dx x = discreteIndex (vcar lo,vcar hi) (vcar dx) (vcar x)
    discreteRange (lo,hi) dx = map (`vcons` vnil) $ discreteRange (vcar lo,vcar hi) (vcar dx)
    discreteRngLength (lo,hi) dx = discreteRngLength (vcar lo,vcar hi) (vcar dx)
    foldRangeM (lo,hi) dx init f = foldRangeM (vcar lo,vcar hi) (vcar dx) init f'
        where f' x acc = f (x `vcons` vnil) acc

instance (Discretizable a, Discretizable (Vec (S n) a)) => Discretizable (Vec (S (S n)) a) where
    discreteIndex (vx, vy) (vd) (vw)
        = discreteIndex (vcar vx,vcar vy) (vcar vd) (vcar vw) +
         (discreteIndex (vcar vx,vcar vy) (vcar vd) (vcar vy) + 1) *
         discreteIndex (vcdr vx,vcdr vy) (vcdr vd) (vcdr vw)
    discreteRange (vx, vy) vd
        = [ vcons xys vxsys | vxsys <- discreteRange (vcdr vx, vcdr vy) (vcdr vd),
                              xys <- discreteRange (vcar vx,vcar vy) (vcar vd)]
    discreteRngLength (vx, vy) (vd) = discreteRngLength (vcar vx,vcar vy) (vcar vd) *
                                      discreteRngLength (vcdr vx,vcdr vy) (vcdr vd)
--    foldRangeM :: Monad m => (a,a) -> a -> b -> (a -> b -> m b) -> m b
    foldRangeM (vlo, vhi) vd  init f =
      let f' x acc =
             foldRangeM (vcdr vlo, vcdr vhi) (vcdr vd) acc $ \ix acc' -> f (x `vcons` ix) acc'
      in foldRangeM (vcar vlo, vcar vhi) (vcar vd) init f'

f v (a, i) = return $ (v:a, i+1)

tf = print =<< foldRangeM (0::Vec Three Int, 2) (1) ([],0) f

tf2 = discreteIndex (0::Vec Three Int, 2) (1) (2 `vcons` 2 `vcons` 2 `vcons` vnil)

{-$ \x-> do
         let thisf ix acc = f (x `vcons` ix) acc
         foldRangeM (vcdr vlo, vcdr vhi) (vcdr vd) init thisf -}

type family Volume a :: *
type instance Volume Time = (Time,Time)
type instance Volume Freq = (Freq,Freq)

--type family FInverses a b :: * -> *

--type instance FInverses Time Freq

--type family FInv a :: *

--type instance FInv Time = Freq
--type instance (FInv t ~ s) => FInv s = t

{-class Inverses a b | a -> b, b-> a where
    invertFwd :: a -> b
    invertBack :: b -> a

instance Inverses Time Freq where
    invertFwd (Time t) = Freq $ recip t
    invertBack (Freq f) = Time $ recip f

instance Inverses Length SpaceFreq where
    invertFwd (Length t) = SpaceFreq $ recip t
    invertBack (SpaceFreq f) = Length $ recip f-}

class HasInverse a where
      type Inv a :: *
      invertValue :: a -> Inv a

instance HasInverse Int where
      type Inv Int = Int
      invertValue = negate -- haha. FIXME

instance HasInverse Time where
      type Inv Time = Freq
      invertValue (Time t) = Freq $ recip t

instance HasInverse Freq where
      type Inv Freq = Time
      invertValue (Freq t) = Time $ recip t

instance HasInverse Length where
      type Inv Length = SpaceFreq
      invertValue (Length t) = SpaceFreq $ recip t

instance HasInverse SpaceFreq where
      type Inv SpaceFreq = Length
      invertValue (SpaceFreq t) = Length $ recip t

instance HasInverse a => HasInverse (Vec n a) where
      type Inv (Vec n a) = Vec n (Inv a)
      invertValue = fmap invertValue

data Signal a b where
    Signal  :: (Discretizable a, {-Ix (Discretized a),-} Storable b, Storable a) =>
               a -> (a,a) -> (Inv a, Inv a) -> SV.Vector b -> Signal a b
    SigFmap :: (b->b') -> Signal a b -> Signal a b'
    SigFun  :: (a->b) -> Signal a b

sigLims :: Signal a b -> (a,a)
sigLims (Signal d1 lims1 invlims1 arr1) = lims1
sigLims (SigFmap _ s) = sigLims s

instance Functor (Signal a) where
    fmap f (SigFmap g s) = SigFmap (f . g) s
    fmap f s = SigFmap f s

sOfVToVofS :: Signal a (Vec n b) -> Vec n (Signal a b)
sOfVToVofS s = undefined

vOfSToSofV :: Vec n (Signal a b) -> Signal a (Vec n b)
vOfSToSofV voss = let sigs = vecToList voss
                  in undefined

type Events a b = [(a,b)]
type Region a b = [(Volume a,b)]

at :: Signal a b -> a -> b
at (Signal delta lims _ arr) x =  arr `SV.index` (discreteIndex lims delta x)
at (SigFmap f s) x = f $ s `at` x
at (SigFun f) x = f x

fill :: (Discretizable a, Storable b, Storable a) =>
        a -> (a,a) -> (Inv a, Inv a) -> (a->b) -> Signal a b
fill delta lims invlims f =
     Signal delta lims invlims $ SV.pack $ map f $ discreteRange lims delta

trans :: Storable c => (a->c) -> Signal a b -> Signal a c
trans f (Signal delta lims ilims arr) = Signal delta lims ilims newarr where
    newarr = SV.pack $ map f $ discreteRange lims delta

trans' :: Storable c => (a->b->c) -> Signal a b -> Signal a c
trans' f (Signal delta lims ilims arr) = Signal delta lims ilims newarr where
    newarr = SV.pack $ map (uncurry f) $ zip (discreteRange lims delta) (SV.unpack arr)
--    newarr = SV.zipWith f (SV.pack $ discreteRange lims delta) (arr)

--DOESNT WORK -- problem may be in foldRangeM implementation
fillIO' :: (Discretizable a, Storable b, Show a, Storable a) =>
          a -> (a,a) -> (Inv a, Inv a) -> (a->IO b) -> IO (Signal a b)
fillIO' delta lims invlims mf = do
  let n = discreteRngLength lims delta
  print n
  stv <- stToIO $ SVST.new_ (n+1800) --FUDGE
  let f ixv nix = do v <- mf ixv
                     --print (ixv, nix)
                     stToIO $ SVST.write stv nix v
                     return $ nix+1
  _ <- foldRangeM lims delta 0 f
  sv <- stToIO $ SVST.unsafeFreeze stv
  return $ Signal delta lims invlims sv

fillIO :: (Discretizable a, Storable b, Show a, Storable a) =>
          a -> (a,a) -> (Inv a, Inv a) -> (a->IO b) -> IO (Signal a b)
fillIO delta lims invlims mf = do
  let rng = discreteRange lims delta
  let n = discreteRngLength lims delta
  xs <- mapM mf rng
  return $ Signal delta lims invlims $ SV.pack xs

zipSignalsWith :: Storable c => (a -> b -> c) -> Signal i a -> Signal i b -> Signal i c
zipSignalsWith f (Signal d1 lims1 invlims1 arr1) (Signal d2 lims2 invlims2 arr2)
               = Signal d1 lims1 invlims2 $ SV.zipWith f arr1 arr2

--toCArray :: Signal a b -> CArray

class X a b where
    fourier :: Signal a b -> Signal (Inv a) b

instance X Time (Complex Double) where
    fourier (Signal d lims invlims arr) = Signal (invertValue d) invlims lims $ SV.pack $ fft $ SV.unpack arr

instance X Freq (Complex Double) where
    fourier (Signal d lims invlims arr) = Signal (invertValue d) invlims lims $ SV.pack $ ifft $ SV.unpack arr

instance X Length (Complex Double) where
    fourier (Signal d lims invlims arr) = Signal (invertValue d) invlims lims $ SV.pack $ fft $ SV.unpack arr

instance X SpaceFreq (Complex Double) where
    fourier (Signal d lims invlims arr) = Signal (invertValue d) invlims lims $ SV.pack $ ifft $ SV.unpack arr

--instance X (Vec Z Length) (Complex Double) where
--    fourier = undefined

instance (Nat n, X (Vec n Length) (Complex Double), Storable (Vec n SpaceFreq),
          Discretizable (Vec n SpaceFreq), Ix (Vec n Int), CA.Shapable (Vec n Int))
    => X (Vec n Length) (Complex Double) where
    fourier s@(Signal d lims invlims arr) =
        let ixdims :: Nat n => Signal (Vec n a) b -> n
            ixdims = undefined
            dims = toInt $ ixdims s
            (oldp, oldn, oldsomething) = SVB.toForeignPtr arr
            cixs = discreteLims s
            carray = unsafePerformIO (CA.unsafeForeignPtrToCArray oldp cixs)
            newcarray = dftN [0..dims] carray
            (newn, newp) = CA.toForeignPtr newcarray
            newarr = SVB.fromForeignPtr newp newn
        in Signal (invertValue d) invlims lims newarr

--instance IsDouble d => Shapable (Vec Z d) where

discreteLims :: (Nat n, IsDouble d) => Signal (Vec n d) a -> (Vec n Int, Vec n Int)
discreteLims s@(Signal d (lo,hi) invlims arr)
    = ( fmap (round . toDouble) $ lo/d,  fmap (round . toDouble) $ hi/d)

instance (Storable (b,c), X a b, X a c) => X a (b,c) where
    fourier s = zipSignalsWith (,) (fourier $ fmap fst s) (fourier $ fmap snd s)

instance (X a b) => X a (Vec n b) where
    fourier sOfVecs =  vOfSToSofV $ fmap fourier $ sOfVToVofS sOfVecs

class XRC a where
    fourierRC :: Signal a Double -> Signal (Inv a) (Complex Double)
    fourierCR :: Signal (Inv a) (Complex Double) -> Signal a Double

max_stars_repo_name: glutamate/space

max_issues_repo_name: glutamate/space

max_forks_repo_name: glutamate/space


text:
{-# LANGUAGE QuantifiedConstraints,
             TupleSections,
             DefaultSignatures,
             FlexibleInstances,
             FlexibleContexts,
             UndecidableInstances,
             ConstraintKinds,
             MultiParamTypeClasses,
             TypeFamilies #-}

module Tensor where
import Data.Complex
import Helper
import EitherTrans
import Control.Monad.Trans.State
import Data.Monoid
import Data.List
import Data.Function
import System.Random
import System.Random.Stateful

-- orphan instance
-- sort complex numbers by their magnitudes
instance RealFloat a => Ord (Complex a) where
    compare a b = compare (magnitude a) (magnitude b)

instance (UniformRange a) => UniformRange (Complex a) where
    uniformRM (lr :+ li, hr :+ hi) g = do
        mr <- uniformRM (lr, hr) g
        mi <- uniformRM (li, hi) g
        return $ mr :+ mi

instance (Uniform a) => Uniform (Complex a) where
    uniformM g = do
        r <- uniformM g
        i <- uniformM g
        return $ r :+ i

{-
    The definition of tensor product should follow conventional tensor product properties.
    Tensor product should form a monoid upto isomorphism
-}

-- higer-kinded definition
-- Polymorphic constraints

class Cmplx c where
    conj :: c -> c
    conj = error "please define conjugate function"

instance (Num a) => Cmplx (Complex a) where
    conj = conjugate

class (Monoid (t a), Cmplx a) => Tensor t a where
    (|*|) :: t a -> t a -> t a
    (|*|) = (<>)
    tconcat :: [t a] -> t a
    tconcat = mconcat
    norm :: Floating a => t a -> t a
    fromList :: Integral b => [a] -> b -> t a
    rank :: Integral b => t a -> b
    isSquare :: t a -> Bool
    conjTr :: t a -> t a
    inner :: t a -> t a -> t a

class (Tensor q a, Num (q a)) => QuantumRegister q a where

    densityOp :: q a
        -> q a

    densityOp x = x * conjTr x

    -- Builds a transformation tensor to apply input gate to selected qubit
    build :: (Integral b, Integral d) => b -- ^ Number of qubits this tensor applies to
        -> [(d, q a)]-- ^ Index / Gate pair (Pre-sorted)
        -> q a -- ^ Default Gate
        -> q a
    build i xs d = tconcat $ go i (sortOnFst xs) d 0
        where
            go i2 [] d2 n = replicate (fromIntegral i2 - n) d2
            go i2 (x:xs) d2 n = replicate (fromIntegral (fst x) - n) d2 ++ (snd x : go i2 xs d2 (fromIntegral $ fst x + 1))

    getProb :: Int -- ^ Index of the qubit
        -> q a -- ^ Qubits
        -> Either String (a, a) -- ^ Resulting Probability
    default getProb :: (Foldable q, Floating a) => Int -> q a -> Either String (a, a) -- ^ Resulting Probability
    getProb i q
        | i < 0 || i >= rank q = Left "Index out of bounds"
        | not (isSquare q) = Left "Invalid Quantum Register"
        | otherwise = do
            a <- applyGate mask0 i q
            b <- applyGate mask1 i q
            return (f a, f b)
        where
            f = getSum . foldMap (Sum . (\x -> conj x * x))

    initQubit :: a -- ^ Constant for |0>
        -> a -- ^ Constant for |1>
        -> Either String (q a) -- ^ Resulting Qubit

    initQubit0 :: q a -- ^ Initialize |0> qubit at z-direction

    initNumQubit0 :: (Integral b) => b -- ^ Number of qubits to initialize
        -> Either String (q a) -- ^ Initialized quantum register

    toQuantumRegister :: [a] -- ^ State vector array to convert
        -> Either String (q a) -- ^ Converted Quantum Register

    collapse :: Int -- ^ Index of the qubit
        -> Bool -- ^ Collapsed State
        -> q a -- ^ Qubits to collapse
        -> Either String (q a) -- ^ Resulting Tensor

    default collapse :: (Floating a) => Int -> Bool -> q a -> Either String (q a)
    collapse i s q
        | i < 0 || i >= rank q = Left "Index out of bounds"
        | not (isSquare q) = Left "Invalid state vector"
        | s = norm <$> applyGate mask1 i q
        | otherwise = norm <$> applyGate mask0 i q

    measure :: (Monad m, StatefulGen g m, UniformRange a) => Int -- ^ Index of the qubit
        -> g -- ^ generator
        -> StateT (q a) (EitherT String m) Bool -- ^ Resulting State Transformer

    default measure :: (Monad m, Num a, Ord a, StatefulGen g m, UniformRange a) => Int -- ^ Index of the qubit
        -> g -- ^ generator
        -> StateT (q a) (EitherT String m) Bool -- ^ Resulting State Transformer
    measure i g = StateT go
        where
            go q
                | i < 0 || i >= rank q = EitherT $ return $ Left "Index out of bounds"
                | not (isSquare q) = EitherT $ return $ Left "Invalid state vector"
                | otherwise = EitherT $ do
                    p <- uniformRM (0, 1) g
                    return $ do
                        t <- getProb i q
                        if p < fst t
                        then (False, ) <$> collapse i False q
                        else (True, ) <$> collapse i True q

    applyGate :: q a -- ^ Gate tensor
        -> Int -- ^ Index
        -> q a -- ^ Quantum Registers
        -> Either String (q a)
    applyGate a i b
        | not (isSquare a) = Left "Invalid Gate Tensor"
        | not (isSquare b) = Left "Invalid Statevector"
        | i < 0 || i >= rank b = Left "Index out of bounds"
        | otherwise = Right $ build (rank b) [(i, a)] pauliId * b

    applyGateAll :: q a -- ^ Gate tensor
        -> q a -- ^ Quantum Registers
        -> Either String (q a)

    applyGateAll a b
        | not (isSquare a) = Left "Invalid Gate Tensor"
        | not (isSquare b) = Left "Invalid Statevector"
        | otherwise = Right $ build (rank b) [] a * b

    applyControl :: q a -- ^ Gate tensor
        -> Int -- ^ Control qubit index
        -> Int -- ^ Apply qubit index
        -> q a -- ^ Quantum Registers
        -> Either String (q a)

    applyControl a ctl i b
        | not (isSquare a) = Left "Invalid gate tensor"
        | not (isSquare b) = Left "Invalid state vector"
        | i < 0 || i >= rank b = Left "Application qubit index out of bounds"
        | ctl < 0 || ctl >= rank b = Left "Control qubit index out of bounds"
        | ctl == i = Left "Control and application qubit cannot be the same"
        | otherwise = Right $ (build (rank b) [(ctl, mask0)] pauliId
            + build (rank b) (sortBy (compare `on` fst) [(ctl, mask1), (i, a)]) pauliId) * b

    illegalPeek :: q a
        -> String

    subSystem :: [Int] -> q a -> Either String (q a)

    isEntangled :: [Int] -> q a -> Either String Bool

    mask0 :: q a
    mask1 :: q a
    pauliId :: q a

class (QuantumRegister q a) => Gates q a where
    pauliX :: q a
    pauliY :: q a
    pauliZ :: q a
    hadamard :: q a

-- Num-only declaration

-- class (forall a. Num a => Monoid (t a)) => Tensor t where
--     (|*|) :: (Num a) => t a -> t a -> t a
--     (|*|) = (<>)
--     norm :: Floating a => t a -> t a
--     fromList :: Integral b => [a] -> b -> t a
--     rank :: Integral b => t a -> b
--     isSquare :: t a -> Bool

-- class (Tensor q, forall a. Num a => Num (q a)) => QuantumRegister q where
--     -- Builds a transformation tensor to apply input gate to selected qubit
--     build :: (RealFloat a, Integral b, Integral c) => b -- ^ Number of qubits this tensor applies to
--         -> [(c, q (Complex a))]-- ^ Index / Gate pair (Pre-sorted)
--         -> q (Complex a) -- ^ Default Gate
--         -> q (Complex a)
--     build i xs d = mconcat $ go i (sortOnFst xs) d 0
--         where
--             go i2 [] d2 n = replicate (fromIntegral i2 - n) d2
--             go i2 (x:xs) d2 n = replicate (fromIntegral (fst x) - n) d2 ++ (snd x : go i2 xs d2 (fromIntegral $ fst x + 1))

--     getProb :: (RealFloat a) => Int -- ^ Index of the qubit
--         -> q (Complex a) -- ^ Qubits
--         -> Either String (a, a) -- ^ Resulting Probability

--     initQubit :: (RealFloat a) => Complex a -- ^ Constant for |0>
--         -> Complex a -- ^ Constant for |1>
--         -> Either String (q (Complex a)) -- ^ Resulting Qubit

--     initQubit0 :: (RealFloat a) => q (Complex a) -- ^ Initialize |0> qubit at z-direction

--     initNumQubit0 :: (RealFloat a, Integral b) => b -- ^ Number of qubits to initialize
--         -> Either String (q (Complex a)) -- ^ Initialized quantum register

--     toQuantumRegister :: (RealFloat a) => [Complex a] -- ^ State vector array to convert
--         -> Either String (q (Complex a)) -- ^ Converted Quantum Register

--     collapse :: (RealFloat a) => Int -- ^ Index of the qubit
--         -> Complex a -- ^ Collapsed State
--         -> q (Complex a) -- ^ Qubits to collapse
--         -> EitherT String IO (q (Complex a)) -- ^ Resulting Tensor

--     measure :: (RealFloat a, Monad m) => Int -- ^ Index of the qubit
--         -> StateT (q (Complex a)) (EitherT String m) (Complex a)
--     -- -> q (Complex a) -- ^ Qubits
--     -- -> Either String ( m (Complex a, q (Complex a))) -- ^ Resulting measurement and the collapsed qubits

--     pauliX :: (RealFloat a) => q (Complex a)
--     pauliY :: (RealFloat a) => q (Complex a)
--     pauliZ :: (RealFloat a) => q (Complex a)
--     pauliId :: (RealFloat a) => q (Complex a)
--     hadamard :: (RealFloat a) => q (Complex a)
--     mask0 :: (RealFloat a) => q (Complex a)
--     mask1 :: (RealFloat a) => q (Complex a)


text:
module Numeric.MixtureModel.Exponential ( -- * General data types
                                          Sample
                                        , SampleIdx
                                        , Samples
                                        , ComponentIdx
                                        , Assignments
                                        , Weight
                                        -- * Exponential parameters
                                        , Rate, Beta
                                        , Exponential(..)
                                        , ComponentParams
                                        , paramFromSamples
                                        , paramsFromAssignments
                                        -- * Exponential distribution
                                        , Prob
                                        , prob
                                        , tauMean, tauVariance
                                        , modelProb
                                        -- * Gibbs sampling
                                        , estimateWeights
                                        , updateAssignments
                                        -- * Score
                                        , scoreAssignments
                                        , maxLikelihoodScore
                                        -- * Classification
                                        , classify
                                        ) where

import           Control.Monad.ST
import           Data.Function (on)
import qualified Data.Vector as VB
import           Data.Vector.Algorithms.Heap
import qualified Data.Vector.Mutable as MV
import qualified Data.Vector.Unboxed as V

import           Numeric.Log hiding (Exp, sum)
import qualified Numeric.Log as Log
import           Numeric.SpecFunctions (logBeta)
import           Statistics.Sample (mean)
import           Numeric.Newton
import           Math.Gamma (gamma)

import           Data.Random hiding (gamma)
import           Data.Random.Distribution.Categorical

type Prob = Log Double
type Sample = Double
type SampleIdx = Int
type ComponentIdx = Int
type Weight = Double

type Rate = Double -- ^ The rate parameter
type Beta = Double -- ^ The stretching parameter
data Exponential = Exp Rate
                 | StretchedExp Rate Beta
                 | FixedExp Rate Beta
                 deriving (Show, Read, Eq)

-- k refers to number of components
-- N refers to number of samples
type Samples = V.Vector Sample                          -- length == N
type Assignments = V.Vector ComponentIdx                -- length == N
type ComponentParams = VB.Vector (Weight, Exponential)  -- length == K

-- | `expProb lambda tau` is the probability of `tau` under Exponential
-- distribution defined by rate `lambda`
prob :: Exponential -> Sample -> Prob
prob _ tau | tau < 0 = error "Exponential distribution undefined for tau<0"
prob (FixedExp lambda 1) tau = prob (Exp lambda) tau
prob (FixedExp lambda beta) tau = prob (StretchedExp lambda beta) tau
prob (StretchedExp lambda 1) tau = prob (Exp lambda) tau
prob (StretchedExp lambda beta) tau =
    Log.Exp $ log beta + (beta-1) * log tau + beta * log lambda - (tau * lambda)**beta
prob (Exp lambda) tau = Log.Exp $ log lambda - lambda * tau

-- | Probability of a sample under a mixture
modelProb :: ComponentParams -> Sample -> Prob
modelProb _ tau | tau < 0 = error "Exponential distribution undefined for tau<0"
modelProb model tau = VB.sum $ VB.map (\(w,e)->prob e tau) model

-- | Mean of the given distribution
tauMean :: Exponential -> Double
tauMean (Exp lambda) = 1 / lambda
tauMean (StretchedExp lambda beta) = gamma (1/beta) / beta / lambda
tauMean (FixedExp lambda beta) = tauMean (StretchedExp lambda beta)

-- | Variance of the given distribution
tauVariance :: Exponential -> Double
tauVariance (Exp lambda) = 1/lambda^2
tauVariance (StretchedExp lambda beta) = 2 * gamma (2/beta) / lambda^2 / beta
tauVariance (FixedExp lambda beta) = tauVariance (StretchedExp lambda beta)

-- | Exponential parameter from samples
paramFromSamples :: Exponential -> V.Vector Sample -> Exponential
paramFromSamples (FixedExp lambda beta) _ = FixedExp lambda beta
paramFromSamples _ v | V.null v = error "Can't estimate parameters without samples"
paramFromSamples (Exp _) v = Exp $ 1 / mean v
paramFromSamples (StretchedExp _ betaOld) v =
    let v' = runST $ do a <- V.thaw v
                        sort a
                        V.freeze a
        n = realToFrac $ V.length v'
        tn = V.head v'
        s beta = V.sum (V.map (\t->t**beta) v') - n*tn**beta
        betaOpt beta = let num = V.sum (V.map (\t->t**beta * log t) v') - n*tn**beta*log tn
                       in num / s beta - V.sum (V.map log v') / n - 1 / beta
        betaOpt' beta = let denom = V.sum (V.map (**beta) v') - n*tn**beta
                            num1 = V.sum (V.map (\t->t**beta * log t) v')
                                   - n*tn**beta * log tn
                            num2 = V.sum (V.map (\t->t**beta * (log t)^2) v')
                                   - n*tn**beta * (log tn)^2
                        in 1/beta^2 - (num1/denom)^2 + num2/denom
        beta' = findRoot 1e-5 betaOpt betaOpt' betaOld
        lambda' = (s beta' / n)**(-1/beta')
    in StretchedExp lambda' beta'

-- | Exponential parameter for component given samples and their
-- component assignments
paramFromAssignments :: Samples -> Assignments -> (ComponentIdx, Exponential) -> Exponential
paramFromAssignments samples assignments (k,p) =
  paramFromSamples p $ V.map snd $ V.filter (\(k',_)->k==k') $ V.zip assignments samples

-- | Exponential parameters for all components given samples and their
-- component assignments
paramsFromAssignments :: Samples -> VB.Vector Exponential -> Assignments -> VB.Vector Exponential
paramsFromAssignments samples params assignments =
  VB.map (paramFromAssignments samples assignments) $ VB.indexed params

-- | Draw a new assignment for a sample given beta parameters
drawAssignment :: ComponentParams -> Sample -> RVar ComponentIdx
drawAssignment params x =
  let probs = map (\(w,p)->realToFrac w * prob p x) $ VB.toList params
  in case filter (isInfinite . ln . fst) $ zip probs [0..] of
        (x:_)     -> return $ snd x
        otherwise -> categorical
                     $ map (\(p,k)->(realToFrac $ p / sum probs :: Double, k))
                     $ zip probs [0..]

-- | `countIndices n v` is the list of counts
countIndices :: Int -> V.Vector Int -> VB.Vector Int
countIndices n v = runST $ do
    accum <- VB.thaw $ VB.replicate n 0
    V.forM_ v $ \k -> do n' <- MV.read accum k
                         MV.write accum k $! n'+1
    VB.freeze accum

-- | Estimate the component weights of a given set of parameters
estimateWeights :: Assignments -> VB.Vector Exponential -> ComponentParams
estimateWeights assignments params =
  let counts = countIndices (VB.length params) assignments
      norm = realToFrac $ V.length assignments
      weights = VB.map (\n->realToFrac n / norm) counts
  in VB.zip weights params

-- | Sample new assignments for observations under given model parameters
updateAssignments :: Samples -> ComponentParams -> RVar Assignments
updateAssignments samples params =
  V.mapM (drawAssignment params) samples

-- | "Likelihood" of sample assignments under given model
-- parameters. Note that the exponential distribution is a density
-- function and as such this will give an unnormalized result unless
-- multiplied by dtau^N
scoreAssignments :: Samples -> ComponentParams -> Assignments -> Prob
scoreAssignments samples params assignments =
    V.product
    $ V.map (\(k,x)->let (w,p) = params VB.! k
                     in realToFrac w * prob p x
            )
    $ V.zip assignments samples

-- | Maximum likelihood classification
classify :: ComponentParams -> Sample -> ComponentIdx
classify params x =
  fst
  $ VB.maximumBy (compare `on` \(_,(w,p))->realToFrac w * prob p x)
  $ VB.indexed params

-- | Score of the maximum likelihood assignment for a set of samples
maxLikelihoodScore :: ComponentParams -> Samples -> Prob
maxLikelihoodScore params samples =
  let assignments = V.map (classify params) samples
  in scoreAssignments samples params assignments

max_stars_repo_name: bgamari/mixture-model

max_issues_repo_name: bgamari/mixture-model

max_forks_repo_name: bgamari/mixture-model


text:
module EchoTest where

import Data.Complex
import Data.Vector.FFT
import qualified Data.Vector.Unboxed as Unboxed
import Echo
import Linear.V3
import Linear.Vector
import Test.Hspec
import Test.QuickCheck

echoTest :: SpecWith ()
echoTest = describe ("Echo signal " ++ (show $ length [0, 1 / fd .. 1e-5])) $ do
  it "no Doppler shift" $ do
    Unboxed.maxIndex (Unboxed.map magnitude harmonicSpectrum) `shouldBe` Unboxed.maxIndex (Unboxed.map magnitude echoSpectrum)
  it "with Doppler shift" $ do
    Unboxed.maxIndex (Unboxed.map magnitude harmonicSpectrum) `shouldNotBe` Unboxed.maxIndex (Unboxed.map magnitude echoSpectrumWithDoppler)

j :: RealFloat a => Complex a
j = 0 :+ 1

f0 = 1000000000.0

fd = 6.325 * f0

harmonicSignal t = exp (- j * (2 * pi * f0 * t :+ 0))

harmonicSpectrum = fft s
  where
    s = Unboxed.fromList $ map harmonicSignal [0, 1 / fd .. 1e-5]

echoSpectrum = fft e
  where
    e = Unboxed.fromList $ map (echo harmonicSignal zero (V3 0 10000 0)) [0, 1 / fd .. 1e-5]

echoSpectrumWithDoppler = fft e
  where
    e = Unboxed.fromList $ map s [0, 1 / fd .. 1e-5]
    s t = echo harmonicSignal (carrier t) (V3 0 10000 0) t

vel :: V3 Double
vel = V3 1000 1000 0

carrier :: Double -> V3 Double
carrier t = zero + vel ^* t

max_stars_repo_name: nihlete/sar-signal

max_issues_repo_name: nihlete/sar-signal

max_forks_repo_name: nihlete/sar-signal


text:
{-# LANGUAGE ScopedTypeVariables #-}

module Clustering.Hierarchical where

import Numeric.LinearAlgebra
import Data.List.Lens
import Control.Lens
import Control.Applicative
import Debug.Trace ( trace )

import Clustering.Util
import Clustering.Data

findClosest2 :: (Ord t) => (a -> a -> t) -> [a] -> Maybe ((Int, Int), t)
findClosest2 _ []  = Nothing
findClosest2 f cls = getidx r1 r2
    where
        r1 = findClosest1s f cls
        r2 = argminval (snd <$>) r1
        getidx :: [Maybe (Int, t)] -> Maybe (Int, Maybe t) -> Maybe ((Int, Int), t)
        getidx l (Just (i, Just v1)) = case l !! i of
          Nothing -> Nothing
          Just (j, v2) -> Just ((i, j), v1)
        getidx l (Just (i, Nothing)) = Nothing
        getidx _ Nothing = Nothing

findClosest1s :: (Ord t) => (a -> a -> t) -> [a] -> [Maybe (Int, t)]
findClosest1s _ [] = []
findClosest1s f ls = reverse $ findClosest1s' [] ls []
    where
        findClosest1s' left [] cl = cl
        findClosest1s' left (x:right) cl =
            findClosest1s' (left ++ [x]) right $
                        selectMMin
                            (findClosest f x left)
                            ((_1 +~ 1 + length left) <$> findClosest f x right)
                            -- shift right index
                        : cl

findClosest :: (Ord t) =>  (a -> a -> t) -> a -> [a] -> Maybe (Int, t)
findClosest _ _ [] = Nothing
findClosest d x ls = argminval (d x) ls

selectMMin :: Ord b => Maybe (a, b) -> Maybe (a, b) -> Maybe (a, b)
selectMMin (Just (n, v1)) (Just (m, v2)) | v2 < v1 = Just (m, v2)
selectMMin (Just (n, v1)) (Just (m, v2)) | v2 > v1 = Just (n, v1)
selectMMin m1 m2 = m1 <|> m2

merge ::(Ord t) => t -> ([a] -> [a] -> t) -> [[a]] -> [[a]]
merge _ _ [] = []
merge m f ls = case findClosest2 f ls of
  Nothing          -> ls
  Just ((i, j), v) -> if v < m
                        then merge m f $ merge2 i j ls
                        else ls

merge2 :: Int -> Int -> [[a]] -> [[a]]
merge2 i j ls = deleteList j $ ix i %~ (\l -> l ++ ls !! j) $ ls

-- | Generalized aggregative scheme using Dendrogram
gasd :: forall a.
    (Dendrogram a -> Dendrogram a -> Float) -- ^ dendrogram distance
    -> (Dendrogram a -> Dendrogram a -> a)  -- ^ representative update
    -> [a] -> Dendrogram a
gasd dd ru xs = agglomerated' 0 $ map dLeaf xs
    where
        agglomerated' :: Int -> [Dendrogram a] -> Dendrogram a
        agglomerated' _ [] = DNil
        agglomerated' _ ds@(d:_) | length ds == 1 = d
        agglomerated' level ds = case findClosest2 dd ds of
            Nothing -> error "Nothing to merge!"
            Just ((i, j), v) -> let (newl, newr) = (ds !! i, ds !!j) in
                agglomerated' (level + 1) $ (DNode (ru newl newr) (level + 1) v newl newr) : dels [i, j] ds

-- | Generalized aggregative scheme
gas ::  Ord t => forall a. ([a] -> [a] -> t) -> [a] -> [[[a]]]
gas cd xs = agglomerate' [map (:[]) xs]
    where
        -- agglomerate' :: [[[a]]] -> [[[a]]]
        agglomerate' [] = []
        agglomerate' [[]] = []
        agglomerate' acs@(c:_) | length c == 1 = acs
        agglomerate' acs@(c:_) = case findClosest2 cd c of
            Nothing -> acs
            Just ((i, j), _) -> agglomerate' $ merge2 i j c:acs

buildSimilarityMtrx :: (a -> a -> Float) -> [a] -> Matrix Float
buildSimilarityMtrx d xs = l >< l $ concatMap (\x -> map (d x) xs) xs
    where
        l = length xs

-- | Matrix update aggregative scheme
muas :: Matrix Float -> (Matrix Float -> Int -> Int -> Int) -> [[Int]]
muas p d = undefined -- TODO optimized agglomerative scheme using matrix update

-- | Modified general divisive scheme
mgds :: (Ord t, Num t, Fractional t) =>
    (a -> a -> t) -- ^ dissimiliarity measure
    -> [a]
    -> [[[a]]]
mgds d [] = []
mgds d xs = divideAll [[xs]]
    where
        divideAll [] = []
        divideAll acs@(c:cs) = if length cs == length xs then acs
                                                         else divideAll $ (c >>= splitCluster) : acs
        splitCluster [] = []
        splitCluster xs = reassign $ outsider xs
                      -- map (((1/fromIntegral (length xs) * ) . sum) . (\x -> map (d x) xs)) xs
        outsider xs = case argmax id $ map (avg . flip map xs . d) xs of
                        Nothing -> error "What happend?!"
                        Just n -> case pop n xs of
                             (as, Nothing) -> ([], as)
                             (as, Just a) -> ([a], as)

        reassign ([], os) = [os]
        reassign (ns, []) = [ns]
        reassign (ns, as) = [ns, as]

max_stars_repo_name: zitkat/clustering

max_issues_repo_name: zitkat/clustering

max_forks_repo_name: zitkat/clustering


text:
{-# LANGUAGE
 TupleSections,
 FlexibleContexts,
 TypeFamilies
 #-}

module DPmixture (
                  stick,
                  dp,
                  dpMixture,
                  dpClusters,
--                  dpMem,
--                  dpMemClusters,
                  posteriorClustersDist
                 ) where

-- Dirichlet Process mixture of Gaussians

import Prelude hiding (sum)

import Data.List hiding (sum)
import Data.Maybe
import Data.Ix (range)
import Numeric.SpecFunctions (logGamma, factorial)
-- import Control.Monad.Memo (memo, startEvalMemoT)
-- import Control.Monad.Trans.Memo.StateCache

import Control.Monad.Bayes.Simple

type NormalInvGamma = (Double,Double,Double,Double)

-- | Prior over cluster parameters used in 'dpMixture'
params :: NormalInvGamma
params = (0,1/10,1,10)
(m,k,a,b) = params

obs :: [Double]
obs = [1.0,1.1,1.2,-1.0,-1.5,-2.0,0.001,0.01,0.005,0.0]

-- | Stick-breaking function.
stick :: MonadDist m => [CustomReal m] -> [a] -> m a
stick (b:breaks) (a:atoms) = do
  stop <- bernoulli b
  if stop then return a else stick breaks atoms

-- | A Dirichlet Process generates a random probability distribution
-- using the stick-breaking representation.
dp :: MonadDist m => CustomReal m -> m a -> m (m a)
dp concentration base = do
  breaks <- sequence $ repeat $ beta 1 concentration
  atoms  <- sequence $ repeat base
  return $ stick breaks atoms

-- | DP mixture example from http://dl.acm.org/citation.cfm?id=2804317
dpMixture :: (MonadBayes m, CustomReal m ~ Double) => m [Int]
dpMixture =
  let
    --lazily generate clusters
    clusters = do
      let atoms = [1..]
      breaks <- sequence $ repeat $ beta 1 1
      let classgen = stick breaks atoms
      vars <- sequence $ repeat $ fmap ((1/) . (*k)) $ gamma a b
      means <- mapM (normal m) vars
      return (classgen, vars, means)
    n = length obs
    --start with no data points
    start = fmap (,[]) clusters

    --add points one by one
    points = foldl build start obs
    build d y = do
      (clusters, rest) <- d
      let (classgen, vars, means) = clusters
      cluster <- classgen
      let mean = means !! cluster
      let var  = vars  !! cluster
      let point = (cluster, var, mean)
      observe (normalDist mean (sqrt var)) y
      return (clusters, point : rest)

  in
   --exctract cluster assignments
   fmap (reverse . map (\(x,_,_) -> x) . snd) points

-- | 'dpMixture' restricted to the number of clusters only
dpClusters :: (MonadBayes m, CustomReal m ~ Double) => m Int
dpClusters = fmap (maximum . normalizePartition) dpMixture

-- | Renames identifiers so that they range from 1 to n
normalizePartition :: [Int] -> [Int]
normalizePartition xs = mapMaybe (`lookup` trans) xs where
  trans = zip unique [1..]
  unique = nub xs

-------------------------------------------------------
-- Strict version using memoization

-- cluster_param :: (MonadDist m, CustomReal m ~ Double) => Int -> m (Double, Double, Double)
-- cluster_param _ = do
--   b <- beta 1 1
--   var <- do
--     t <- gamma a b
--     return (1 / (k*t))
--   mean <- normal m var
--   return (b, mean, var)
--
-- dpMem :: (MonadBayes m, CustomReal m ~ Double) => m [Int]
-- dpMem = startEvalMemoT $ mapM process_point obs where
--   process_point x = do
--     (c, mean, var) <- get_cluster
--     observe (Continuous (Normal mean (sqrt var))) x
--     return c
--   get_cluster = stick 0 where
--     stick c = do
--       (b, mean, var) <- memo cluster_param c
--       stop <- bernoulli b
--       if stop then return (c, mean, var) else stick (c+1)
--
-- dpMemClusters :: (MonadBayes m, CustomReal m ~ Double) => m Int
-- dpMemClusters = fmap (maximum . normalizePartition) dpMem

----------------------------------------------------------
-- Exact posterior

-- | Posterior for parameters of a cluster
posterior :: NormalInvGamma -> [Double] -> NormalInvGamma
posterior (m,k,a,b) xs = (m',k',a',b') where
    m' = (k * m + n * x) / k'
    k' = k + n
    a' = a + n / 2
    b' = b + (m * m * k + x2 - m' * m' * k') / 2 --s / 2 + (k * n * (x - m) ^ 2) / (2 * (k + n))
    n  = fromIntegral $ length xs
    x  = sum xs / n  -- mean xs
    x2 = sum $ zipWith (*) xs xs  -- sum xs^2
    s  = sum $ zipWith (*) d d where  -- pvar xs
            d = map (+ (- x)) xs

-- | Model evidence for a cluster
evidence :: NormalInvGamma -> [Double] -> Double
evidence prior xs = exp (logGamma a' - logGamma a) * (b ** a / b' ** a') * sqrt (k / k') * c ^ n where
    (m,k,a,b)     = prior
    (m',k',a',b') = posterior prior xs
    c             = 1 / (2 * (sqrt pi))
    n             = length xs

-- | CRP prior over partitions
crp :: [Int] -> Double
crp labels = (product $ map factor counts) / factorial n where
    factor x = factorial (x - 1)
    n = length labels
    counts = map (\x -> length (filter (== x) labels)) clusters
    clusters = nub labels

-- | DP mixture likelihood of a partition
mixture :: NormalInvGamma -> [Double] -> [Int] -> Double
mixture prior xs ks = (product $ map (evidence prior) clusters) * crp ks where
  n = length xs
  clusters = map extract labels
  labels = nub ks
  extract n = map fst $ filter ((== n) . snd) points
  points = zip xs ks

-- | Exact posterior for a partition
exact :: NormalInvGamma -> [Double] -> [Int] -> Double
exact prior xs = let
  z = sum $ map (mixture prior xs) (partitions n)
  n = length xs
  f partition = mixture prior xs partition / z
  in
   f

-- | Exact posterior over the number of clusters
exactClusters :: NormalInvGamma -> [Double] -> Int -> Double
exactClusters prior xs = let
  partitionScore = exact prior xs
  results = map score [1..n]
  score k = sum $ map partitionScore $ filter ((== k) . maximum) $ partitions n
  n = length xs
  lookup i = results !! (i - 1)
  in
   lookup

-- | All possible partitions
partitions :: Int -> [[Int]]
partitions n = generate n 1 where
  generate 0 _ = [[]]
  generate n k = do
    x <- range (1,k)          -- range is inclusive on both ends
    let k' = max k (x+1)      -- next empty cluster
    xs <- generate (n-1) k'
    return (x:xs)

-- | Posterior over the number of clusters
posteriorClustersDist :: (MonadDist m, CustomReal m ~ Double) => m Int
posteriorClustersDist = categorical $ zip ns $ map lookup ns where
    ns = [1 .. length obs]
    lookup = exactClusters params obs

--------------------------------------
-- Instances for memoization
-- type instance CustomReal (StateCache c m) = CustomReal m
--
-- instance MonadDist m => MonadDist (StateCache c m) where
--   primitive = lift . primitive
--
-- instance MonadBayes m => MonadBayes (StateCache c m) where
--   factor = lift . factor

max_stars_repo_name: eligottlieb/monad-bayes

max_issues_repo_name: eligottlieb/monad-bayes

max_forks_repo_name: eligottlieb/monad-bayes


text:
module Main where

import Net
import Data.Traversable
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Devel
import qualified Data.Vector.Storable as V
import System.Random
import Data.List
import Control.Monad

testNet :: DoubleNet -> [Vector Double] -> [Vector Double] -> Double
testNet n ins ans =
  let outs :: [Vector Double]
      outs = map (runNet n) ins
  in sum (zipWith (\x y -> norm_1 $ add x (scale (-1) y)) outs ans)

rankNets :: [DoubleNet] -> IO [(Double,DoubleNet)]
rankNets nets =
  do values <- mapM (\_ -> do x <- randomRIO (0,10)
                              return $ vector [x]) [1..100]
     let answers = map sin values
         scoredNets = map (\x -> (testNet x values answers,x)) nets
     return $ sortOn fst scoredNets

plotSome :: DoubleNet -> IO ()
plotSome n =
  do let pts = map (\x -> (x/1000)*10) $ [1..1000]
     mapM_ (\x -> putStrLn $ (show $ runNet n (vector [x])) ++ " " ++ (show $ (sin x)/2)) pts

sinFunc = DiffFunction sin (\x y -> y*cos x)

errFunc :: Double -> DifferentiableFunction (Vector Double) Double
errFunc c = DiffFunction (\v -> ((atIndex v 0) - c)^2) (\v w -> 2*((atIndex v 0) - c)*(atIndex w 0))

iterateNet :: Int -> DoubleNet -> IO DoubleNet
iterateNet 0 net = return net
iterateNet steps net =
  do let n = 10
     xs <- replicateM n $ randomRIO (0,20)

     --let x = 0.5
     let loop loopSteps currNet =
           if loopSteps ==0 then
             return currNet
           else
             do grads <- forM xs (\x -> do let (val,err,grad) = runAndDiffWithError currNet (vector [x]) (errFunc $ (sin x)/2)
                                         --print (x,val,(sin x)/2,err)
                                           return grad)
                let grad = foldl1 (zipWith (\a b -> (add b a))) grads
                let newNet = updateNet net grad (-0.0005)
                loop (loopSteps - 1) newNet
     loop 5 net >>= iterateNet (steps-1)

main :: IO ()
main = do randNet <- randomSimpleNet 1 1 100
          finalNet <- iterateNet 10000 randNet
          plotSome finalNet

max_stars_repo_name: jkelleyy/sinenet

max_issues_repo_name: jkelleyy/sinenet

max_forks_repo_name: jkelleyy/sinenet


text:
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE FlexibleContexts #-}

module Neural.Matrix(
    initNet, updateNet, feedLayer, forwardPass, propagate,
    backPropagate, relu, softmax, categoricalCrossEntropy
) where

import Control.Monad (zipWithM)
import Data.List (foldl', scanl', transpose)
import Control.Arrow ((&&&))
import Numeric.AD (grad)
import Numeric.AD.Mode (auto)
import Numeric.AD.Internal.Reverse (Reverse, Tape)
import Data.Reflection (Reifies)
import Numeric.LinearAlgebra

type Biases = Vector R
type Weights = Matrix R
type Activation = [R] -> [R]
type Layer = (Activation, Biases, Weights)

parameters :: Layer -> (Biases, Weights)
parameters (_, bs, ws) = (bs, ws)

weights :: Layer -> Weights
weights (_, _, ws) = ws

initNet :: R -> [Int] -> [Activation] -> IO [Layer]
initNet b szs@(_:lns) acts = layer <$> zipWithM randn szs lns
  where layer = zip3 acts $ vector . flip replicate b <$> lns

relu :: Activation
relu = map $ max 0

relu' :: Activation
relu' = map f
  where f x | x < 0     = 0
            | otherwise = 1

categoricalCrossEntropy :: RealFloat a => [a] -> [a] -> a
categoricalCrossEntropy ys = negate . sum . zipWith (*) ys . map log . softmax

softmax :: RealFloat a => [a] -> [a]
softmax ts = let exps = exp <$> ts in map (/ sum exps) exps

feedLayer :: (Vector R, Vector R) -> Layer -> (Vector R, Vector R)
feedLayer (_, an) (aF, bs, ws) = (id &&& fromList . aF . toList) $ bs + an <# ws

forwardPass :: Vector R -> [Layer] -> [(Vector R, Vector R)]
forwardPass = scanl' feedLayer . (id &&& id)

propagate :: (Vector R, Matrix R) -> Vector R -> Vector R
propagate (zn, wn) dldn = wn #> dldn * dadz
  where dadz = fromList . relu' . toList $ zn

backPropagate :: Vector R -> [Vector R] -> [Matrix R] -> [Vector R]
backPropagate dldL tls lws = tail $ scanr propagate dldL $ zip tls lws

updateLayer :: R -> (Vector R, Vector R, Layer) -> Layer
updateLayer s (a, e, (aF, bs, ws)) = (aF, newBs, newWs)
  where newWs = ws - scalar s * outer a e
        newBs = bs - scalar s * e

updateNet :: (forall a. Reifies a Tape => [Reverse a R]
                                       -> [Reverse a R]
                                       -> Reverse a R) ->
             R -> Vector R -> Vector R -> [Layer] -> [Layer]
updateNet c s xs ys ls = map (updateLayer s) $ zip3 (init ans) dls ls
  where fps  = forwardPass xs ls
        ans  = snd <$> fps
        yhat = auto <$> (toList . last) ans
        dldL = fromList $ grad (c (auto <$> (toList ys))) yhat
        dls  = backPropagate dldL (fst <$> init fps) (weights <$> ls)

max_stars_repo_name: Ultramann/hnet

max_issues_repo_name: Ultramann/hnet

max_forks_repo_name: Ultramann/hnet


text:
-- |
-- Copyright: (C) 2013 Amgen, Inc.
--

{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeFamilies #-}

module Data.Vector.SEXP.Base where

import Control.Memory.Region

import Foreign.R.Type
import Foreign.R (SEXP, SomeSEXP)

import Data.Singletons (SingI)

import Data.Complex (Complex)
import Data.Word (Word8)
import Data.Int (Int32)
import Foreign.Storable (Storable)

-- | Function from R types to the types of the representations of each element
-- in the vector.
type family ElemRep s (a :: SEXPTYPE) where
  ElemRep s 'Char    = Word8
  ElemRep s 'Logical = Logical
  ElemRep s 'Int     = Int32
  ElemRep s 'Real    = Double
  ElemRep s 'Complex = Complex Double
  ElemRep s 'String  = SEXP s 'Char
  ElemRep s 'Vector  = SomeSEXP s
  ElemRep s 'Expr    = SomeSEXP s
  ElemRep s 'Raw     = Word8

-- | 'ElemRep' in the form of a relation, for convenience.
type E s a b = ElemRep s a ~ b

-- | Constraint synonym for all operations on vectors.
type VECTOR s ty a = (Storable a, IsVector ty, SingI ty)

-- | Constraint synonym for all operations on vectors.
type SVECTOR ty a = (Storable a, IsVector ty, SingI ty, ElemRep V ty ~ a)

max_stars_repo_name: reactormonk/HaskellR

max_issues_repo_name: reactormonk/HaskellR

max_forks_repo_name: reactormonk/HaskellR