text:
module Lib where

import qualified Numeric.LinearAlgebra as Matrix
import qualified Numeric.LinearAlgebra.Data as HM

import Control.Monad.Random
import Control.Monad.Random.Class

import Control.DeepSeq

import Control.Parallel.Strategies (parListChunk, rpar, rdeepseq, rseq, runEval, runEvalIO)

import Data.List (foldl', sort, sortBy, maximumBy, minimumBy)

import GHC.Conc (numCapabilities)

import qualified Data.ByteString.Lazy as BS (writeFile, readFile)
import Data.Serialize (encodeLazy, decodeLazy)

data Player = X | O deriving (Show, Ord, Eq, Enum, Read)

type Grid   = (Player, [Maybe Player])
data Result = Undecided | Tie | Win Player deriving Show

data Tree a = Node a [Tree a] deriving (Show)

showGrid :: Grid -> String
showGrid (_, xs) = " " ++ x00 ++ " | " ++ x01 ++ " | " ++ x02 ++ "\n" ++
                   "-----------" ++ "\n" ++
                   " " ++ x10 ++ " | " ++ x11 ++ " | " ++ x12 ++ "\n" ++
                   "-----------" ++ "\n" ++
                   " " ++ x20 ++ " | " ++ x21 ++ " | " ++ x22 ++ "\n"
    where
        [x00, x01, x02, x10, x11, x12, x20, x21, x22] =
            map (\x -> case x of
                          Nothing -> " "
                          Just X  -> "X"
                          Just O  -> "O") xs

moves :: (Grid, Result) -> [(Grid, Result)]
moves (g@(_, xs), Undecided) = gs
    where
        l  = length $ filter id $ map isNothing xs
        gs = map (\i -> move i g) [0..l-1]
moves _ = []

gameTree :: (Grid, Result) -> Tree (Grid, Result)
gameTree n = Node n fs
    where
        fs = map gameTree $ moves n

maxMin :: RandomGen g => Tree (Grid, Result) -> Rand g (Maybe ((Grid, Result), Int))
maxMin (Node _ []) = return Nothing
maxMin (Node n fs) = case mms of
                        [] -> return Nothing
                        _  -> do
                               i <- getRandomR (0, length maxScores - 1)
                               return (Just $ maxScores !! i)
    where
        mms = map minMax' fs
        maxScore  = maximum $ map snd mms
        maxScores = filter ((==maxScore) . snd) mms

maxMin' :: Tree (Grid, Result) -> ((Grid, Result), Int)
maxMin' (Node n@(_, Tie) []) = (n, 0)
maxMin' (Node n@(_, Win O) []) = (n, -1)
maxMin' (Node n@(_, Win X) []) = (n, 1)
maxMin' (Node (_, _) []) = error "Undecided Leaf"
maxMin' (Node n fs) = (n, snd $ maximumBy (\(_, r1) (_, r2) -> r1 `compare` r2) $ map minMax' fs)

minMax :: RandomGen g => Tree (Grid, Result) -> Rand g (Maybe ((Grid, Result), Int))
minMax (Node _ []) = return Nothing
minMax (Node n fs) = case mms of
                        [] -> return Nothing
                        _  -> do
                               i <- getRandomR (0, length minScores - 1)
                               return (Just $ minScores !! i)
    where
        mms = map maxMin' fs
        minScore  = minimum $ map snd mms
        minScores = filter ((==minScore) . snd) mms

minMax' :: Tree (Grid, Result) -> ((Grid, Result), Int)
minMax' (Node n@(_, Tie) []) = (n, 0)
minMax' (Node n@(_, Win O) []) = (n, -1)
minMax' (Node n@(_, Win X) []) = (n, 1)
minMax' (Node (_, _) []) = error "Undecided Leaf"
minMax' (Node n fs) = (n, snd $ minimumBy (\(_, r1) (_, r2) -> r1 `compare` r2) $ map maxMin' fs)

moveMM :: RandomGen g => (Grid, Result) -> Rand g (Maybe (Grid, Result))
moveMM gr@((X, _), _) = fmap (fmap fst) . maxMin $ gameTree gr
moveMM gr@((O, _), _) = fmap (fmap fst) . minMax $ gameTree gr

empty :: Grid
empty = (X, take 9 $ repeat Nothing)

chunksOf :: Int -> [a] -> [[a]]
chunksOf n [] = []
chunksOf n xs = take n xs : chunksOf n (drop n xs)

rows :: [a] -> [[a]]
rows xs      = rs ++
               map (\i -> map (!!i) rs) [0,1,2] ++
               pure (zipWith (!!) rs [0,1,2]) ++
               pure (zipWith (!!) rs [2,1,0])
    where
        rs = chunksOf 3 xs

allSame :: Eq a => [a] -> (Bool, Maybe a)
allSame []     = (True, Nothing)
allSame (x:xs) = case all (==x) xs of
                   True  -> (True, Just x)
                   False -> (False, Nothing)

isNothing :: Maybe a -> Bool
isNothing Nothing = True
isNothing _       = False

gridResult :: Grid -> Result
gridResult (_, xs)
        | not (null winners) = Win $ head winners
        | undecided          = Undecided
        | otherwise          = Tie
    where
        rs = rows xs
        rresults = map rowResult rs
        winners = map (\(Just x) -> x) $
                  filter (not . isNothing) $
                  map (\res -> case res of
                                 Win p -> Just p
                                 _     -> Nothing)
                  rresults
        undecided = any (\res -> case res of
                                 Undecided -> True
                                 _         -> False)
                  rresults

rowResult :: [Maybe Player] -> Result
rowResult r = case any isNothing r of
                True  -> Undecided
                False -> case allSame r of
                           (True, Just (Just p)) -> Win p
                           _                     -> Tie

emptyFields :: [Maybe a] -> [Int]
emptyFields xs = map snd $
                 filter (isNothing . fst) $
                 zip xs [0..]

move :: Int -> Grid -> (Grid, Result)
move m (p, xs) = ((nextPlayer, nextGrid), res)
    where
        nextPlayer = case p of
                       X -> O
                       O -> X
        empties = emptyFields xs
        numEmpties = length empties
        m' = m `mod` numEmpties
        pos = empties !! m'
        nextGrid = take pos xs ++ [Just p] ++ drop (pos+1) xs
        res = gridResult (nextPlayer, nextGrid)

moveAI :: DenseData -> Grid -> (Grid, Result)
moveAI net g = move m g
    where
        l = fromGrid g
        m = predict (net, [relu, softmax]) l

battle :: DenseData -> DenseData -> Maybe Player
battle pX pO = battle' pX pO empty

battle' :: DenseData -> DenseData -> Grid -> Maybe Player
battle' pX pO g@(p, _) = case moveAI ai g of
                            (_, Tie)        -> Nothing
                            (_, Win w)      -> Just w
                            (g', Undecided) -> battle' pX pO g'
    where
        ai = case p of
                X -> pX
                O -> pO

battleMM' :: RandomGen g => (Grid, Result) -> Rand g [(Grid, Result)]
battleMM' g = do
                g' <- moveMM g
                case g' of
                  Nothing  -> return [g]
                  Just g'' -> fmap (g:) $ battleMM' g''

battleMM :: RandomGen g => Rand g [(Grid, Result)]
battleMM = battleMM' (empty, Undecided)

loses :: DenseData -> DenseData -> Bool
loses p1 p2 = case w1 of
                Just O -> True
                _      -> case w2 of
                            Just X -> True
                            _      -> False
    where
        w1 = battle p1 p2
        w2 = battle p2 p1

beats :: DenseData -> DenseData -> Bool
beats p1 p2 = case w1 of
                Just X -> case w2 of
                            Nothing -> True
                            Just O  -> True
                            _       -> False
                _      -> False
    where
        w1 = battle p1 p2
        w2 = battle p2 p1

type Layer  = HM.Vector HM.R
type Bias   = HM.Vector HM.R
type Weight = HM.Matrix HM.R
type Activation = Layer -> Layer

type DenseData = [(Bias, Weight)]
type DenseNet  = (DenseData, [Activation])

type CrossingData = [(Bias, Bias, Weight, Weight)]

toFile :: String -> DenseData -> IO ()
toFile fileName net = BS.writeFile fileName $
                      encodeLazy $
                      fmap (\(b, w) -> (HM.toList b, HM.toLists w)) net

fromFile :: String -> IO DenseData
fromFile fileName = fmap (fmap (\(b, w) -> (HM.fromList b, HM.fromLists w)) . (\(Right x) -> x) . decodeLazy) $
                    BS.readFile fileName

randBias :: RandomGen g => Int -> Rand g Bias
randBias n = (pure . HM.fromList . take n) =<< getRandomRs (-1, 1)

randWeight :: RandomGen g => Int -> Int -> Rand g Weight
randWeight m n = (pure . (m HM.>< n)) =<< getRandomRs (-1, 1)

repeatM :: Monad m => Int -> m a -> m [a]
repeatM 0 m = pure []
repeatM n m =
    do
        x  <- m
        xs <- repeatM (n-1) m
        pure (x:xs)

toBinary :: (Num a, Ord a) => a -> a
toBinary x
    | x < 0 = 0
    | otherwise = 1

invert :: Num a => a -> a
invert x = 1-x

randCrossing :: RandomGen g => Rand g CrossingData
randCrossing =
  do
    ~[(b1, w01), (b2, w12)] <- randPlayer
    let b1' = HM.cmap toBinary b1
    let b2' = HM.cmap toBinary b2
    let w01' = HM.cmap toBinary w01
    let w12' = HM.cmap toBinary w12
    return [(b1', HM.cmap invert b1', w01', HM.cmap invert w01'),
            (b2', HM.cmap invert b2', w12', HM.cmap invert w12')]

randPlayer :: RandomGen g => Rand g DenseData
randPlayer =
  do
    b1 <- randBias 18
    b2 <- randBias 9
    w01 <- randWeight 18 18
    w12 <- randWeight 9 18
    return [(b1, w01), (b2, w12)]

fromGrid :: Grid -> Layer
fromGrid (_, xs) = HM.fromList $ lX ++ lO
    where
        lX = map (\x -> case x of
                          Just X -> 1
                          _      -> 0)
             xs
        lO = map (\x -> case x of
                          Just O -> 1
                          _      -> 0)
             xs

dense :: Bias -> Weight -> Activation -> Layer -> Layer
dense b w a l = a $ b + (w Matrix.#> l)

relu :: Activation
relu = HM.cmap relu'
    where
        relu' x = if x < 0 then 0 else x

softmax :: Activation
softmax l = HM.cmap (/s) exps
    where
        m    = HM.maxElement l
        exps = HM.cmap (\x -> exp (x - m)) l
        s    = sum $ HM.toList exps

runDenseNet :: DenseNet -> Layer -> Layer
runDenseNet (dat, acts) l =
    foldl' (\l' ((b, w), a) -> dense b w a l') l $ zip dat acts

predict :: DenseNet -> Layer -> Int
predict net = HM.maxIndex . runDenseNet net

select :: [DenseData] -> (Int, [DenseData])
select ps = (maximum scores,
             take n' $
             map fst $
             sortBy (\(_, score1) (_, score2) -> score2 `compare` score1) $
             zip ps scores)
    where
        scores = runEval $
                 parListChunk (n `div` numCapabilities) rdeepseq $
                 fmap (\x -> length $ filter not $ fmap (loses x) ps) ps
        n = length ps
        n' = n `div` 100

breed2' :: DenseData -> DenseData -> CrossingData -> DenseData
breed2' m f cd = map (\((b, w), (b', w'), (cb, cb', cw, cw')) -> (b*cb + b' * cb', w*cw + w'*cw')) $ zip3 m f cd

breed2 :: RandomGen g => DenseData -> DenseData -> Rand g [DenseData]
breed2 m f =
    do
        cds <- repeatM 10 randCrossing
        return $ fmap (breed2' m f) cds

breed :: RandomGen g => [DenseData] -> Rand g [DenseData]
breed ps = fmap concat $ sequence $ breed2 <$> ps <*> ps

initPopulation :: RandomGen g => Int -> Rand g [DenseData]
initPopulation n = repeatM n randPlayer

someFunc :: IO ()
someFunc = do
            ps <- evalRandIO $ initPopulation 1000
            go ps
    where
        go ps = let (score, fittest) = select ps
                in do
                     putStrLn $ show score ++ "\t(" ++ show (length ps) ++ ")"
                     if score == 1000
                       then return () >> (toFile "weights.dat" $ head fittest)
                       else go =<< evalRandIO (breed fittest)

max_stars_repo_name: nilsalex/genetic-ttt

max_issues_repo_name: nilsalex/genetic-ttt

max_forks_repo_name: nilsalex/genetic-ttt


text:
{-# LANGUAGE ForeignFunctionInterface #-}
{-|
Module      : Grenade.Layers.Internal.Transpose
Description : Functions to quickly transpose a matrix
Maintainer  : Theo Charalambous
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Internal.Transpose (
  transpose4d
) where

import qualified Data.Vector.Storable        as U (unsafeFromForeignPtr0,
                                                   unsafeToForeignPtr0)

import           Foreign                     (mallocForeignPtrArray,
                                              withForeignPtr)
import           Foreign.C.Types             (CInt)
import           Foreign.Ptr                 (Ptr)
import           Numeric.LinearAlgebra       (Matrix, Vector, cmap, flatten)
import qualified Numeric.LinearAlgebra.Data  as D
import qualified Numeric.LinearAlgebra.Devel as U
import           System.IO.Unsafe            (unsafePerformIO)

import           Grenade.Types

-- | Transpose a 4d matrix, similarly to the tranpose function in numpy
transpose4d :: [Int]           -- ^ original dimensions
            -> Vector RealNum  -- ^ permutation vector
            -> Matrix RealNum  -- ^ input matrix
            -> Matrix RealNum  -- ^ transposed matrix
transpose4d dims@[n, c, h, w] permsV m
  = let outMatSize = n * c * h * w
        outW       = dims !! (round $ permsV D.! 3)
        vec        = flatten m
        permsV'    = cmap (U.fi . round) permsV  :: Vector CInt
        dimsV      = D.fromList $ map U.fi dims  :: Vector CInt
    in unsafePerformIO $ do
      outPtr        <- mallocForeignPtrArray outMatSize
      let (inPtr, _) = U.unsafeToForeignPtr0 vec
          (pPtr, _)  = U.unsafeToForeignPtr0 permsV'
          (dPtr, _)  = U.unsafeToForeignPtr0 dimsV

      withForeignPtr inPtr $ \inPtr' ->
        withForeignPtr pPtr $ \pPtr' ->
          withForeignPtr dPtr $ \dPtr' ->
            withForeignPtr outPtr $ \outPtr' ->
              transpose_4d inPtr' dPtr' pPtr' outPtr'

      let matVec = U.unsafeFromForeignPtr0 outPtr outMatSize
      return (U.matrixFromVector U.RowMajor (div outMatSize outW) outW matVec)
{-# INLINE transpose4d #-}

foreign import ccall unsafe
    transpose_4d
      :: Ptr RealNum -> Ptr CInt -> Ptr CInt -> Ptr RealNum -> IO ()

ext: hs

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE DataKinds                  #-}
{-# LANGUAGE TypeOperators              #-}

module DFT
  ( testSuite
  ) where

import Data.Complex
import Data.Mod.Word
import Data.Poly.Semiring (UPoly, unPoly, toPoly, dft, inverseDft, dftMult)
import qualified Data.Vector.Unboxed as U
import GHC.TypeNats (KnownNat, natVal, type (+), type (^))
import Test.Tasty
import Test.Tasty.QuickCheck hiding (scale, numTests)

import Dense ()

testSuite :: TestTree
testSuite = testGroup "DFT"
  [ testGroup "dft matches reference"
    [ dftMatchesRef (0 :: Mod (2 ^ 0 + 1))
    , dftMatchesRef (2 :: Mod (2 ^ 1 + 1))
    , dftMatchesRef (2 :: Mod (2 ^ 2 + 1))
    , dftMatchesRef (3 :: Mod (2 ^ 4 + 1))
    , dftMatchesRef (3 :: Mod (2 ^ 8 + 1))
    ]
  , testGroup "dft is invertible"
    [ dftIsInvertible (0 :: Mod (2 ^ 0 + 1))
    , dftIsInvertible (2 :: Mod (2 ^ 1 + 1))
    , dftIsInvertible (2 :: Mod (2 ^ 2 + 1))
    , dftIsInvertible (3 :: Mod (2 ^ 4 + 1))
    , dftIsInvertible (3 :: Mod (2 ^ 8 + 1))
    ]
  , testProperty "dftMult matches reference" dftMultMatchesRef
  ]

dftMatchesRef :: KnownNat n1 => Mod n1 -> TestTree
dftMatchesRef primRoot = testProperty (show n) $ do
  xs <- U.replicateM n arbitrary
  pure $ dft primRoot xs === dftRef primRoot xs
  where
    n = fromIntegral (natVal primRoot - 1)

dftRef :: (Num a, U.Unbox a) => a -> U.Vector a -> U.Vector a
dftRef primRoot xs = U.generate (U.length xs) $
  \k -> sum (map (\j -> xs U.! j * primRoot ^ (j * k)) [0..n-1])
  where
    n = U.length xs

dftIsInvertible :: KnownNat n1 => Mod n1 -> TestTree
dftIsInvertible primRoot = testProperty (show n) $ do
  xs <- U.replicateM n arbitrary
  let ys = dft primRoot xs
      zs = inverseDft primRoot ys
  pure $ xs === zs
  where
    n = fromIntegral (natVal primRoot - 1)

dftMultMatchesRef :: UPoly Int -> UPoly Int -> Property
dftMultMatchesRef xs ys = zs === dftZs
  where
    xs', ys', dftZs' :: UPoly (Complex Double)
    xs' = toPoly $ U.map fromIntegral $ unPoly xs
    ys' = toPoly $ U.map fromIntegral $ unPoly ys
    dftZs' = dftMult (\k -> cis (2 * pi / fromIntegral k)) xs' ys'

    zs, dftZs :: UPoly (Complex Int)
    zs  = toPoly $ U.map (:+ 0) $ unPoly $ xs * ys
    dftZs  = toPoly $ U.map (\(x :+ y) -> round x :+ round y) $ unPoly dftZs'

max_stars_repo_name: Bodigrim/poly

max_issues_repo_name: Bodigrim/poly

max_forks_repo_name: Bodigrim/poly


text:
{- |
Raw, unnormalized multi-dimensional versions of the transforms in @fftw@.

Note that the forwards and backwards transforms of this module are not actually
inverses.  For example, @run idft (run dft v) /= v@ in general.

For more information on the individual transforms, see
<http://www.fftw.org/fftw3_doc/What-FFTW-Really-Computes.html>.

@since 0.2
-}

module Numeric.FFT.Vector.Unnormalized.Multi
  (
        -- * Creating and executing 'Plan's
        run,
        plan,
        execute,
        -- * Complex-to-complex transforms
        dft,
        idft,
        -- * Real-to-complex transforms
        dftR2C,
        dftC2R,
  ) where

import Numeric.FFT.Vector.Base
import Foreign
import Foreign.C
import Data.Complex

#include <fftw3.h>

-- | Whether the complex fft is forwards or backwards.
type CDirection = CInt

foreign import ccall unsafe fftw_plan_dft
    :: CInt -> Ptr CInt -> Ptr (Complex Double) -> Ptr (Complex Double)
        -> CDirection -> CFlags -> IO (Ptr CPlan)

foreign import ccall unsafe fftw_plan_dft_r2c
    :: CInt -> Ptr CInt -> Ptr Double -> Ptr (Complex Double) -> CFlags
        -> IO (Ptr CPlan)

foreign import ccall unsafe fftw_plan_dft_c2r
    :: CInt -> Ptr CInt -> Ptr (Complex Double) -> Ptr Double -> CFlags
        -> IO (Ptr CPlan)

dftND :: CDirection -> TransformND (Complex Double) (Complex Double)
dftND d = TransformND
  { inputSizeND = id
  , outputSizeND = id
  , creationSizeFromInputND = id
  , makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft rk dims a b d
  , normalizationND = const id
  }

-- | A forward discrete Fourier transform.  The output and input sizes are the same (@n@).
dft :: TransformND (Complex Double) (Complex Double)
dft = dftND (#const FFTW_FORWARD)

-- | A backward discrete Fourier transform.  The output and input sizes are the same (@n@).
idft :: TransformND (Complex Double) (Complex Double)
idft = dftND (#const FFTW_BACKWARD)

-- | A forward discrete Fourier transform with real data.  If the input size is @n0 * ... * nk@,
-- the output size will be @n0 * ... * nk \`div\` 2 + 1@.
dftR2C :: TransformND Double (Complex Double)
dftR2C = TransformND {
              inputSizeND = id,
              outputSizeND = \n -> n `div` 2 + 1,
              creationSizeFromInputND = id,
              makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft_r2c rk dims a b,
              normalizationND = const id
          }

-- | A backward discrete Fourier transform which produces real data.
--
-- This 'Transform' behaves differently than the others:
--
--  - Calling @plan dftC2R n@ creates a 'Plan' whose /output/ size is @n@, and whose
--    /input/ size is @n \`div\` 2 + 1@.
--
--  - If @length v == n@, then @length (run dftC2R v) == 2*(n-1)@.
dftC2R :: TransformND (Complex Double) Double
dftC2R = TransformND {
            inputSizeND = \n -> n `div` 2 + 1,
            outputSizeND = id,
            creationSizeFromInputND = \n -> 2 * (n-1),
            makePlanND = \rk dims a b -> withPlanner . fftw_plan_dft_c2r rk dims a b,
            normalizationND = const id
        }

ext: hsc

max_stars_repo_name: TravisWhitaker/vector-fftw

max_issues_repo_name: TravisWhitaker/vector-fftw

max_forks_repo_name: TravisWhitaker/vector-fftw


text:
--{-# LANGUAGE UndecidableInstances,
--             FlexibleInstances,
--             FlexibleContexts,
--             TypeFamilies,
--             ScopedTypeVariables #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Numeric.Signal.Noise
-- Copyright   :  (c) Alexander Vivian Hugh McPhail 2010, 2014, 2015
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  uses Concurrency
--
-- Noise generation functions
--
-----------------------------------------------------------------------------

--             IncoherentInstances,

module Numeric.Signal.Noise (
                             pinkNoise
                            , spatialNoise
                            , powerNoise
                            ) where

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

--import qualified Numeric.Signal as S

--import Complex

--import qualified Data.Array.IArray as I
--import Data.Ix

--import Data.Word

--import System.IO.Unsafe(unsafePerformIO)

--import qualified Data.List as L

--import Data.Binary

--import Foreign.Storable

import Numeric.LinearAlgebra

import qualified Data.Vector.Generic as GV

import qualified Numeric.GSL.Fourier as F

--import qualified Numeric.GSL.Histogram as H
--import qualified Numeric.GSL.Histogram2D as H2

--import qualified Numeric.Statistics.Information as SI

import Prelude hiding(filter)

--import Control.Monad(replicateM)

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

-- | The method is briefly descirbed in Lennon, J.L. "Red-shifts and red
-- herrings in geographical ecology", Ecography, Vol. 23, p101-113 (2000)
--
-- Matlab version Written by Jon Yearsley  1 May 2004
-- j.yearsley@macaulay.ac.uk
--
-- Creates 1/f scale invariant spatial noise
spatialNoise :: Double    -- ^  β: spectral distribution
                         --    0: White noise
                         --   -1: Pink noise
                         --   -2: Brownian noise
             -> Int -> Int -- ^ matrix dimensions
             -> Int       -- ^ random seed
             -> Matrix Double
spatialNoise b r' c' s = let c = fromIntegral c'
                             r = fromIntegral r'
                             pre_x = linspace c' (0::Double,c-1)
                             post_x = linspace c' (c,1)
                             freq_x = GV.map (/c) $ vjoin [pre_x,post_x]
                             u = fromRows (replicate (2*r') freq_x)
                             pre_y = linspace r' (0::Double,r-1)
                             post_y = linspace r' (r,1)
                             freq_y = GV.map (/c) $ vjoin [pre_y,post_y]
                             v = fromColumns (replicate (2*c') freq_y)
                             s_f = cmap (**(b/2)) ((u**2) + (v**2))
                             s_f' = cmap (\x -> if isInfinite x then 0 else x) s_f
                             phi = reshape (2*c') (randomVector s Uniform (4*r'*c'))
                         in subMatrix (1,1) (r',c') $ fst $ fromComplex $ fromRows $ map F.ifft $ toRows $ ((complex $ s_f'**0.5) * (toComplex (cos(2*pi*phi),sin(2*pi*phi))))

-- | 1/f scale invariant noise
pinkNoise ::
    Double    -- ^  β: spectral distribution
              --    0: White noise
              --   -1: Pink noise
              --   -2: Brownian (red) noise
    -> Int       -- ^ samples
    -> Int       -- ^ random seed
    -> Vector Double
pinkNoise b s r = let pre = linspace s (0::Double,fromIntegral (s-1))
                      post = linspace s (fromIntegral s,1)
                      freq = vjoin [pre/(fromIntegral s),post/(fromIntegral s)]
                      s_f = GV.map (**(b/2)) (freq**2)
                      s_f' = GV.map (\x -> if isInfinite x then 0 else x) s_f
                      phi = randomVector r Uniform (2*s)
                  in subVector 0 s $ fst $ fromComplex $ F.ifft ((complex $ s_f'**0.5) * (toComplex (cos(2*pi*phi),sin(2*pi*phi))))

-- | generate noise from a power spectrum
powerNoise :: Vector Double   -- ^ the power spectrum
           -> Int             -- ^ random seed
           -> Vector Double
powerNoise psd r = let ln = GV.length psd
                       freq = vjoin [fromList [0],psd, (fromList . reverse . tail . toList) psd]
                       phi = randomVector r Uniform (2*ln)
                   in (fromIntegral ln) * (subVector 0 (ln-1) $ fst $ fromComplex $ F.ifft ((complex $ freq) * (toComplex (cos(2*pi*phi),sin(2*pi*phi)))))

max_stars_repo_name: amcphail/hsignal

max_issues_repo_name: amcphail/hsignal

max_forks_repo_name: amcphail/hsignal


text:
--Haskell Quil compiler
--Copyright Laurence Emms 2018

module Lib (compile) where

import Data.Complex
import Data.List
import Register
import Instruction
import ClassicalCircuit
import QuantumCircuit
import Util

testRXCircuit :: (Floating a, Show a, Ord a) => Circuit a
testRXCircuit = let parameters@[Left a] = [Left (MetaQubitRegister "a")]
                    instructions = [(RX (ComplexConstant (5.0 :+ 10.0)) a)]
                in Circuit "TESTRX" parameters instructions

hadamardRotateXYZ :: (Floating a, Show a, Ord a) => Complex a -> Circuit a
hadamardRotateXYZ rotation
    = let parameters@[Left a, Left b, Right r, Left z] = [Left (MetaQubitRegister "a"), Left (MetaQubitRegister "b"), Right (MetaRegister "r"), Left (MetaQubitRegister "z")]
          instructions = [(Hadamard a),
                          (MeasureOut a r)] ++
                          (ifC "HROTXYZTHEN0" "HROTXYZEND0"
                               r
                               [(RX (ComplexConstant rotation) z)]
                               ([(Hadamard b),
                                 (MeasureOut b r)] ++
                                (ifC "HROTXYZTHEN1" "HROTXYZEND1"
                                     r
                                     [(RY (ComplexConstant rotation) z)]
                                     [(RZ (ComplexConstant rotation) z)])))
      in Circuit "HROTXYZ" parameters instructions

compile :: IO ()
compile = putStrLn "Compiling quantum executable" >>
          --Test registers
          putStrLn (show $ QubitRegister 10) >>
          putStrLn (show $ Register 5) >>
          putStrLn (show $ ComplexConstant (6.0 :+ 7.0)) >>
          putStrLn (show $ ComplexConstant ((-8.0) :+ 1.0)) >>
          putStrLn (show $ ComplexConstant (2.0 :+ (-3.0))) >>
          putStrLn (show $ ComplexConstant ((-4.0) :+ (-6.0))) >>
          --Test instructions
          putStrLn (show $ CNot (QubitRegister 0) (QubitRegister 1)) >>
          putStrLn (show $ PSwap (ComplexConstant (5.0 :+ (-3.2))) (QubitRegister 0) (QubitRegister 1)) >>
          putStrLn (show $ Measure (QubitRegister 4)) >>
          putStrLn (show $ MeasureOut (QubitRegister 4) (Register 5)) >>
          --Test circuits
          putStrLn (show $ DefCircuit bell) >>
          putStrLn (show $ CallCircuit bell [Left (QubitRegister 5), Left (QubitRegister 3)]) >>
          putStrLn (show $ DefCircuit testRXCircuit) >>
          putStrLn (show $ CallCircuit testRXCircuit [Left (QubitRegister 1)]) >>
          putStrLn (show $ DefCircuit xor) >>
          putStrLn (show $ CallCircuit xor [Right (Register 0), Right (Register 1), Right (Register 2)]) >>
          putStrLn (show $ DefCircuit halfAdder) >>
          putStrLn (show $ CallCircuit halfAdder [Right (Register 0), Right (Register 1), Right (Register 2), Right (Register 3)]) >>
          putStrLn (show $ DefCircuit adder) >>
          putStrLn (show $ CallCircuit adder [Right (Register 0), Right (Register 1), Right (Register 2), Right (Register 3), Right (Register 4)]) >>
          --Load the integer 53 into registers [0-31]
          putStrLn (intercalate "\n" (fmap show (loadIntToRegisters 53 0))) >>
          --Load the integer 18 into registers [32-63]
          putStrLn (intercalate "\n" (fmap show (loadIntToRegisters 18 32))) >>
          --Add [0-31] + [32-63] into [64-95]
          let carry = 128
              temp = 129
          in putStrLn (intercalate "\n" (fmap show (addRegisters carry temp 0 32 64))) >>
          putStrLn (show $ DefCircuit (hadamardRotateXYZ (0.70710678118 :+ 0.70710678118))) >>
          putStrLn (show $ CallCircuit (hadamardRotateXYZ (0.70710678118 :+ 0.70710678118)) [Left (QubitRegister 0), Left (QubitRegister 1), Right (Register 0), Left (QubitRegister 2)])

max_stars_repo_name: WhatTheFunctional/Hasquil

max_issues_repo_name: WhatTheFunctional/Hasquil

max_forks_repo_name: WhatTheFunctional/Hasquil


text:
module CommonFunctions where

import Data.Complex

-- f(z) = z^2 + c
f1 :: Complex Double -> Complex Double -> Complex Double
f1 c z = z**2 + c

-- below is valid haskell
-- f1 :: Double -> Double -> Double
-- f1 a b = a**2 + b


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

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

This module contains functions for displaying the
state of a spin-1/2 particle or other quantum two-level
system as a point on the Bloch Sphere.
-}

module Physics.Learn.BlochSphere
    ( VisObj
    , toPos
    , ketToPos
    , staticBlochSphere
    , displayStaticState
    , animatedBlochSphere
    , simulateBlochSphere
    , simulateBlochSphereK
    , stateProp
    , statePropK
    , evolutionBlochSphere
    , evolutionBlochSphereK
    , hamRabi
    )
    where

import qualified Physics.Learn.QuantumMat as M
import qualified Physics.Learn.Ket as K
import Physics.Learn.Ket
    ( Ket
    , Operator
    , (<>)
    , dagger
    )
import Numeric.LinearAlgebra
    ( Vector
    , Matrix
    , C
    , iC
--    , (<>)  -- matrix multiplication
--    , (|>)  -- vector definition
    , (!)   -- vector element access
    , (><)  -- matrix definition
    , scale
    , size
    )
import Data.Complex
    ( Complex(..)
    , conjugate
    , realPart
    , imagPart
    )
import Physics.Learn
    ( Position
    , v3FromPos
    , cart
    )
import SpatialMath
    ( Euler(..)
    )
import Vis
    ( VisObject(..)
    , Flavour(..)
    , Options(..)
    , Camera0(..)
    , defaultOpts
    , display
    , simulate
    , blue
    , red
    )
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

{-
3 ways to specify the state of a spin-1/2 particle:
Vector C
Ket
Position  (Bloch vector)

2 ways to specify a Hamiltonian:
Matrix C
Operator

3 choices for Vis' world:
(Float, Vector C)
(Float, Ket)
(Float, Position)
-}

-- | A Vis object.
type VisObj = VisObject Double

-- | Convert a 2x1 complex state vector for a qubit
--   into Bloch (x,y,z) coordinates.
toPos :: Vector C -> Position
toPos v
    = if size v /= 2
      then error "toPos only for size 2 vectors"
      else let z1 = v ! 0
               z2 = v ! 1
           in cart (2 * realPart (conjugate z1 * z2))
                   (2 * imagPart (conjugate z1 * z2))
                   (realPart (conjugate z1 * z1 - conjugate z2 * z2))

-- | Convert a qubit ket
--   into Bloch (x,y,z) coordinates.
ketToPos :: Ket -> Position
ketToPos psi
    = if K.dim psi /= 2
      then error "ketToPos only for qubit kets"
      else let z1 = dagger K.zp <> psi
               z2 = dagger K.zm <> psi
           in cart (2 * realPart (conjugate z1 * z2))
                   (2 * imagPart (conjugate z1 * z2))
                   (realPart (conjugate z1 * z1 - conjugate z2 * z2))

-- | A static 'VisObj' for the state of a qubit.
staticBlochSphere :: Position -> VisObj
staticBlochSphere r
    = RotEulerDeg (Euler 270 0 0) $ RotEulerDeg (Euler 0 180 0) $
      VisObjects [ Sphere 1 Wireframe blue
                 , Trans (v3FromPos r) (Sphere 0.05 Solid red)
                 ]

displayStaticBlochSphere :: Position -> IO ()
displayStaticBlochSphere r
    = display myOptions (staticBlochSphere r)

-- | Display a qubit state vector as a point on the Bloch Sphere.
displayStaticState :: Vector C -> IO ()
displayStaticState = displayStaticBlochSphere . toPos

-- | Given a Bloch vector as a function of time,
--   return a 'VisObj' as a function of time.
animatedBlochSphere :: (Double -> Position) -> (Float -> VisObj)
animatedBlochSphere f
    = staticBlochSphere . f . realToFrac

-- | Given a sample rate, initial qubit state vector, and
--   state propagation function, produce a simulation.
--   The 'Float' in the state propagation function is the time
--   since the beginning of the simulation.
simulateBlochSphere :: Double -> Vector C -> (Float -> (Float,Vector C) -> (Float,Vector C)) -> IO ()
simulateBlochSphere sampleRate initial statePropFunc
    = simulate myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc

-- | Given a sample rate, initial qubit state ket, and
--   state propagation function, produce a simulation.
--   The 'Float' in the state propagation function is the time
--   since the beginning of the simulation.
simulateBlochSphereK :: Double -> Ket -> (Float -> (Float,Ket) -> (Float,Ket)) -> IO ()
simulateBlochSphereK sampleRate initial statePropFuncK
    = simulate myOptions sampleRate (0,initial) (staticBlochSphere . ketToPos . snd) statePropFuncK

{-
-- | Given a sample rate, initial qubit state vector, and
--   state propagation function, produce a simulation.
--   The 'Float' in the state propagation function is the time
--   since the beginning of the simulation.
playBlochSphere :: Double -> Vector C -> (Float -> (Float,Vector C) -> (Float,Vector C)) -> IO ()
playBlochSphere sampleRate initial statePropFunc
    = play myOptions sampleRate (0,initial) (staticBlochSphere . toPos . snd) statePropFunc
-}

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

-- | Produce a state propagation function from a time-dependent Hamiltonian.
statePropK :: (Double -> Operator) -> Float -> (Float,Ket) -> (Float,Ket)
statePropK ham tNew (tOld,psi)
    = (tNew, K.timeEv (realToFrac dt) (ham tMid) psi)
      where
        dt = tNew - tOld
        tMid = realToFrac $ (tNew + tOld) / 2

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

-- | Given an initial qubit ket and a time-dependent Hamiltonian,
--   produce a visualization.
evolutionBlochSphereK :: Ket -> (Double -> Operator) -> IO ()
evolutionBlochSphereK psi0 ham
    = simulateBlochSphereK 0.01 psi0 (statePropK ham)

myOptions :: Options
myOptions = defaultOpts {optWindowName = "Bloch Sphere"
                        ,optInitialCamera = Just (Camera0 75 20 4)}

{-
staticBz1 :: IO ()
staticBz1 = evolutionBlochSphere M.xp (const (scale 0.9 M.sz))

staticBz2 :: IO ()
staticBz2 = evolutionBlochSphere ((2|>) [(cos (pi / 8)), (sin (pi / 8))]) (const M.sz)

staticBy1 :: IO ()
staticBy1 = evolutionBlochSphere M.xp (const M.sy)
-}

-- | Hamiltonian for nuclear magnetic resonance.
--   Explain omega0, omegaR, omega.
hamRabi :: Double ->  Double ->  Double ->  Double -> Matrix C
hamRabi omega0 omegaR omega t
    = let h11 = omega0 :+ 0
          h12 = (omegaR :+ 0) * exp (-iC * ((omega * t) :+ 0))
      in scale (1/2) $ (2><2) [h11, h12, (conjugate h12), (-h11)]

-- need to scale time

-- a pi pulse

max_stars_repo_name: walck/learn-physics

max_issues_repo_name: walck/learn-physics

max_forks_repo_name: walck/learn-physics


text:
module Impulse.Impulse where

import           Control.Monad             as M
import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Binary
import           Data.Complex
import           Data.List                 as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           STC.CompletionField
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types
import           Utils.Array

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:scale0FreqsStr:scaleFreqsStr:initDistStr:histFilePath:alphaStr:pinwheelFlagStr:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      initDist = read initDistStr :: [R2S1RPPoint]
      scale0Freq = read scale0FreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      scale0Freqs = [-scale0Freq .. scale0Freq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      alpha = read alphaStr :: Double
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numThread = read numThreadStr :: Int
      sourceDist = L.take 1 initDist
      sinkDist = L.drop 1 initDist
      folderPath = "output/test/Impulse"
      (R2S1RPPoint (_, _, _, s)) = L.head sourceDist
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          theta0Freqs
          thetaFreqs
          scale0Freqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scale0Freqs
             , L.length theta0Freqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      theta0Freqs
      scale0Freqs
  createDirectoryIfMissing True folderPath
  let arr4d =
        R.slice
          arrR2Z2T0S0
          (Z :. All :. All :. (L.length theta0Freqs - 1) :.
           (L.length scale0Freqs - 1) :.
           All :.
           All)
  createDirectoryIfMissing True (folderPath </> "GreensR2Z2T0S0")
  -- MP.mapM_
  --   (\(i, j) ->
  --      plotImageRepaComplex
  --        (folderPath </> "GreensR2Z2T0S0" </> show (i + 1) L.++ "_" L.++
  --         show (j + 1) L.++
  --         ".png") .
  --      ImageRepa 8 .
  --      computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.slice arr4d $
  --      (Z :. i :. j :. All :. All))
  --   [ (i, j)
  --   | i <- [0 .. (L.length thetaFreqs) - 1]
  --   , j <- [0 .. (L.length scaleFreqs) - 1]
  --   ]
  plan <- makeR2Z2T0S0Plan emptyPlan arrR2Z2T0S0
  sourceDistArr <-
    computeInitialDistributionR2T0S0
      plan
      numPoint
      numPoint
      theta0Freqs
      scale0Freqs
      maxScale
      -- initDist
      sourceDist
  arrR2Z2T0S0F <- dftR2Z2T0S0 plan . computeS . makeFilter2D $ arrR2Z2T0S0
  -- Source field
  sourceArr <- convolveR2T0S0 plan arrR2Z2T0S0F sourceDistArr
  sourceR2Z2 <- R.sumP . R.sumS . rotateR2Z2T0S0Array $ sourceArr
  sourceField <-
    fmap (computeS . R.extend (Z :. (1 :: Int) :. All :. All)) .
    R.sumP .
    R.sumS .
    rotate4D .
    rotate4D . r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
    sourceR2Z2
  plotImageRepaComplex (folderPath </> printf "Source%.0f.png" s) . ImageRepa 8 $
    sourceField

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE RecordWildCards #-}
{-
The shape of a triangle can be represented by two explicit points (ie vectors)
and a third implicit point at (0,0). A triangle defined by three points can be
converted by subtracting one of the points from the others:
  shape (p1,p2,p3) = (p2-p1, p3-p1)
This keeps the shape intact but drops the position of the triangle.

If P and Q are two such triangles, there exists a matrix A that transforms
P into Q: AP = Q
The 'A' matrix can be calculated by multiplying Q and the inverse of P:
  A = QP_inv

So, if we have a set of triangles, we can drop their positional information
and compute their shape transformation matrices. But can we do it in reverse?
Can we start with a set of transformation matrices (which don't give us the
positions of the triangles) and compute a set of points that represent triangles
with the right transformation?
In other words, if we have a set of triangles (with absolute position) and
a transformation matrix for each triangle, can we compute new absolute positions
such that AP=Q?

Solving this is not straightforward as moving a single point will change the
transformation of multiple triangles. Fortunately, though, this problem can be
expression as linear equations and solved quickly with LAPACK.

To express the problem as a linear system, we need 4 equations for each
triangle and 2 unknowns for each point.

The system has this shape: Mx=B
We're trying to find 'x' which are the positions of the new points.
'M' is a matrix that transforms points into transformation matrices.
'B' is the target transformation matrices.

With a bit of rewriting, we can reduce 'A = QP_inv' to 4 linear equations
with 6 unknowns. The 6 unknowns are the x and y positions for the three corners
in a triangle.
This is just another way of calculating the transformation matrix 'A' and we
can verify it by looking at the output of:
  * computeA src dst
  * applyCoeff (coeffOfB src) dst
The output should be the same because they're both computing 'A'.

'applyCoeff' simply multiplies the coefficients by a 6x1 matrix containing
the x and y positions for 'dst'.

So, this means we have roughly the right shape: Mx=B. M is the coefficients
we found for 'src', 'x' is dst, and 'B' is the transformation matrix that
turns 'src' into 'dst'.
We can get the coefficients for each triangle and place them together in a
single large matrix. Then we can put the desired transformation matrices in 'B'
and ask a linear solver to find 'x' which satisfies the equations.

Phew, step one is done. We now have the machinery for turning arbitrary
transformation matrices into proper, connected triangles. Next step is figuring
out which rotational matrices to use. At the two extremes, we have the identity
matrices (representing our starting shape) and the matrices computed from the
target triangles (represetning out target shape). How can these two matrices be
interpolated for a smooth morph?

There are several ways of interpolating matrices. Linear interpolation is
possible. But, it looks prettier if we separate rotation from shear. Then we
can use spherical linear interpolation on the rotation and linear interpolation
on the shear.

-}
module Reanimate.Morph.Rigid where

import           Data.Foldable (toList)
import           Data.Vector                   (Vector)
import qualified Data.Vector                   as V
import           Linear.Quaternion
import           Linear.Vector
import           Linear.V2
import           Linear.V3
import qualified Numeric.LinearAlgebra         as Matrix
import           Numeric.LinearAlgebra.HMatrix (GMatrix, Matrix,
                                                toLists, (!), (><))

type P = V2 Double
type Trig = (P,P,P)
type RelTrig = (Int,Int,Int)

data Mesh = Mesh
  { meshPointsA :: Vector P
  , meshPointsB :: Vector P
  , meshOutline :: Vector Int
  -- , meshSteiner :: Vector Int
  , meshTriangles :: Vector RelTrig }

-- applyA (computeA a b) a = b + some_constant_translation
-- A = Q P_inv
computeA :: Trig -> Trig -> Matrix Double
computeA p q = matQ <> Matrix.inv matP
  where
    matP = trigToMatrix p
    matQ = trigToMatrix q

-- A = UDV = U(VV_t)DV = (UV)(V_tDV) = RS
-- R = UV
-- S = V_tDV
computeRS :: Matrix Double -> (Matrix Double, Matrix Double)
computeRS a = (r, s)
  where
    (u,d,vt) = Matrix.svd a
    v = Matrix.tr vt
    r = u <> v
    s = vt <> Matrix.diag d <> v

matrixToQuaternion :: Matrix Double -> Quaternion Double
matrixToQuaternion r =  q
  where
    w = 0.5 * sqrt (1 + r!0!0 + r!1!1 + 1)
    z = 1/(4*w) * (r!0!1 - r!1!0)
    q = Quaternion w (V3 0 0 z)

quaternionToMatrix :: Quaternion Double -> Matrix Double
quaternionToMatrix q = (2><2)
    [ 1 - 2*(qj*qj + qk*qk), 2*(qi*qj+qk*qr)
    , 2*(qi*qj - qk*qr), 1 - 2*(qi*qi + qk*qk) ]
  where
    [qr,qi,qj,qk] = toList q

-- A = R((1-t)I + tS)
computeA_RSt :: Matrix Double -> Matrix Double -> Double -> Matrix Double
computeA_RSt r s t = r_t <> (realToFrac (1-t) * Matrix.ident 2 + realToFrac t * s)
  where
    i = Quaternion 1 0
    q = slerp i (matrixToQuaternion r) t
    r_t = quaternionToMatrix q

applyA :: Matrix Double -> Trig -> Trig
applyA a p =
    case toLists (a <> matP) of
      [ [x1, x2], [y1, y2] ] -> (V2 0 0, V2 x1 y1, V2 x2 y2)
      _ -> error "invalid matrix"
  where
    matP = trigToMatrix p

coeffOfB :: Trig -> Matrix Double
coeffOfB p = (4><6)
    [ -a0-a2,      0, a0,   0, a2,  0
    , -a1-a3,      0, a1,   0, a3,  0
    ,      0, -a0-a2,  0,  a0,  0, a2
    ,      0, -a1-a3,  0,  a1,  0, a3 ]
  where
    [[a0,a1],[a2,a3]] = toLists (Matrix.inv (trigToMatrix p))

applyCoeff :: Matrix Double -> Trig -> Matrix Double
applyCoeff b (V2 x1 y1, V2 x2 y2, V2 x3 y3) = b <> matQ
  where
    matQ = (6><1) [ x1, y1, x2, y2, x3, y3]

trigToMatrix :: Trig -> Matrix Double
trigToMatrix (p1,p2,p3) = matP
  where
    V2 p12x p12y = p2-p1
    V2 p13x p13y = p3-p1
    matP = (2><2)
      [p12x, p13x
      ,p12y, p13y]

data Prep = Prep
  { prepPivot    :: (P, P)
  , prepPointsA  :: Vector P
  , prepPointsB  :: Vector P
  , prepRS       :: Vector (Matrix Double, Matrix Double)
  , prepRSRev    :: Vector (Matrix Double, Matrix Double)
  , prepUToB     :: GMatrix
  }

symmetric :: Bool
symmetric = True

prepare :: Mesh -> Prep
prepare Mesh{..} = Prep
    { prepPivot = (aOrigin, bOrigin)
    , prepPointsA =
        V.map (subtract aOrigin) $
        V.take pivotIdx meshPointsA <> V.drop (pivotIdx+1) meshPointsA
    , prepPointsB =
        V.map (subtract bOrigin) $
        V.take pivotIdx meshPointsB <> V.drop (pivotIdx+1) meshPointsB
    , prepRS = rsList
    , prepRSRev = rsRevList
    , prepUToB = Matrix.mkSparse uToB }
  where
    aOrigin = meshPointsA V.! pivotIdx
    bOrigin = meshPointsB V.! pivotIdx
    pivotIdx = case V.head meshTriangles of
        (a,_,_) -> a
    mkAbs p (a,b,c) = (p V.! a,p V.! b,p V.! c)
    absATrigs = V.map (mkAbs meshPointsA) meshTriangles
    absBTrigs = V.map (mkAbs meshPointsB) meshTriangles
    aList = V.zipWith computeA absATrigs absBTrigs
    rsList = V.map computeRS aList
    revList = V.zipWith computeA absBTrigs absATrigs
    rsRevList = V.map computeRS revList
    n = length meshTriangles
    -- nT = if symmetric then n*2 else n
    -- pivotToB = ((nT*4)><2) $ concat
    --   [ [ fromMaybe 0 (lookup (x,pivotIdx*2) bigM)
    --     , fromMaybe 0 (lookup (x,pivotIdx*2+1) bigM)]
    --   | x <- [0..(nT*4)-1] ]
    uToB =
      [ ((x,y-2),key) | ((x,y),key) <- bigM, y /= pivotIdx*2 && y /= pivotIdx*2+1 ]
    bigM =
      concat (zipWith worker [0..] (V.toList meshTriangles)) ++
      if symmetric
        then concat $ zipWith workerRev [n..] (V.toList meshTriangles)
        else []
    worker i src@(a,b,c) = concat $
      let effs = coeffOfB (mkAbs meshPointsA src) in
      [ [((i*4+h, e*2), effs!h!(j*2))
        ,((i*4+h, e*2+1), effs!h!(j*2+1))]
      | h <- [0..3]
      , (e,j) <- zip [a,b,c] [0..]
      ]
    workerRev i dst@(a,b,c) = concat $
      let effs = coeffOfB (mkAbs meshPointsB dst) in
      [ [((i*4+h, e*2), effs!h!(j*2))
        ,((i*4+h, e*2+1), effs!h!(j*2+1))]
      | h <- [0..3]
      , (e,j) <- zip [a,b,c] [0..]
      ]

interpolate :: Prep -> Double -> Vector P
interpolate Prep{..} t = V.fromList $
    pivot : worker (Matrix.toList solution)
  where
    -- solution = Matrix.cgSolve False prepUToB b
    solution = Matrix.cgx $ last solutions
    solutions = Matrix.cgSolve'
      False
      1e-9
      1e-9
      1000
      prepUToB
      b
      (Matrix.fromList $ concat [ [x,y] | V2 x y <- V.toList target ])
      -- 0
    target = if t < 0.5
      then prepPointsA
      else prepPointsB
    worker (x:y:xs) = V2 x y ^+^ pivot : worker xs
    worker _ = []
    pivot = case prepPivot of
      (src, dst) -> lerp t dst src
    n = V.length prepRS
    b = Matrix.vector $
      [ concat (toLists a)!!j
      | i <- [0..n-1]
      , let (r,s) = prepRS V.! i
      , let a = computeA_RSt r s t
      , j <- [0..3]
      ] ++
      [ concat (toLists a)!!j
      | symmetric
      , i <- [0..n-1]
      , let (r,s) = prepRSRev V.! i
      , let a = computeA_RSt r s (1-t)
      , j <- [0..3]
      ]

max_stars_repo_name: TristanCacqueray/reanimate

max_issues_repo_name: TristanCacqueray/reanimate

max_forks_repo_name: TristanCacqueray/reanimate


text:
module Parser
  ( readExpr
  ) where

import           Control.Monad.Except
import           Data.Char                     (toLower)
import           Data.Complex                  (Complex (..), realPart)
import           Data.Ratio
import           Eval
import           Numeric                       (readHex, readOct)
import           Prelude                       hiding (toRational)
import           Syntax
import           Text.ParserCombinators.Parsec hiding (spaces)

-- parse toma: Un Parser como primer parametro
--             Un String como nombre de la entrada, para los errores.
--             La entrada
readExpr :: String -> ThrowsError LispVal
readExpr input = case parse parseExpr "lisp" input of
                   Left err -> throwError $ Parser err
                   Right value  -> return value

parseExpr :: Parser LispVal
parseExpr = try parseNumber
            <|> try parseString
            <|> try  parseBool
            <|> try parseCharacter
            <|> try parseAtom
            <|> try parseQuoted
            <|> try parseListLike

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

parseBool :: Parser LispVal
parseBool = do _ <- char '#'
               c <- oneOf "tf"
               (return . Bool) $ case c of
                                   't' -> True
                                   'f' -> False
-- Parseo de strings
parseString :: Parser LispVal
parseString = do _ <- char '"'
                 x <- many (parseEscapedChar <|> noneOf ['\\', '"'])
                 _ <- char '"'
                 return (String x)

parseEscapedChar :: Parser Char
parseEscapedChar = do _ <- char '\\'
                      c <- oneOf ['\\', '"', 'n', 'r', 't']
                      return $ case c of
                                 '\\' -> c
                                 '"'  -> c
                                 'n'  -> '\n'
                                 'r'  -> '\r'
                                 't'  -> '\t'

-- Parseo de caracteres
parseCharacter :: Parser LispVal
parseCharacter = do _ <- string "#\\"
                    s <- many letter
                    (return . Character) $ case map toLower s of
                                             "spaces"  -> ' '
                                             ""        -> ' '
                                             "newline" -> '\n'
                                             [x]       -> x

-- Parseo de tipos numericos
parseNumber :: Parser LispVal
parseNumber = try parseComplexNumber
           <|> try parseFloat
           <|> try parseRationalNumber
           <|> try parseInteger

parseInteger :: Parser LispVal
parseInteger = try parseIntWithPrefix <|> try parsePlainInt

parseIntWithPrefix :: Parser LispVal
parseIntWithPrefix = do _ <- char '#'
                        try parseBin <|> try parseDec <|> try parseOct <|> try parseHex

parsePlainInt :: Parser LispVal
parsePlainInt = (Number . read) <$> many1 digit

parseFloat :: Parser LispVal
parseFloat = do int <- many1 digit
                _       <- char '.'
                decimal <- many1 digit
                let number = int ++ "." ++ decimal
                return $ (Float . read) number

-- Toma un numero de lisp y devuelve un numero de haskell
-- asumo que ese numero puede es real (un double)
fromNumber :: LispVal -> Double
fromNumber (Number x)  = fromIntegral x
fromNumber (Float x)    = x
fromNumber (Rational x) = fromRational x
fromNumber (Complex x)  = realPart x

parseComplexNumber :: Parser LispVal
parseComplexNumber = do real <- fromNumber <$> (try parseFloat <|> try parseRationalNumber <|> try parseInteger)
                        _ <- oneOf "+-"
                        imag <- fromNumber <$> (try parseFloat <|> try parseRationalNumber <|> try parseInteger)
                        return $ Complex (real :+ imag)

parseRationalNumber :: Parser LispVal
parseRationalNumber = do p <- read <$> many1 digit
                         _ <- char '/'
                         q <- read <$> many1 digit
                         return $ Rational (p % q)

-- Parseo de numeros en diferentes bases numericas
readBin :: Char -> Int
readBin '0' = 0
readBin '1' = 1

fromBinary :: String -> Int
fromBinary s = sum $ map (\(i,x) -> i*(2^x)) $ zip [0..] $ map readBin (reverse s)

parseBin :: Parser LispVal
parseBin = do _ <- char 'b'
              number <- many1 (oneOf "01")
              return $ (Number . fromBinary) number

parseDec :: Parser LispVal
parseDec = do _ <- char 'd'
              number <- many1 digit
              return $ (Number . read) number

parseOct :: Parser LispVal
parseOct = do _ <- char 'o'
              number <- many1 (oneOf "01234567")
              return $ (Number . fst . head . readOct) number

parseHex :: Parser LispVal
parseHex = do _ <- char 'x'
              number <- many1 (oneOf "0123456789abcdefABCDEF")
              return $ (Number . fst . head . readHex) number

-- Parseo de listas y listas puntuadas
parseListLike :: Parser LispVal
parseListLike = do _ <- char '('
                   x <- try parseList <|> parseDottedList
                   _ <- char ')'
                   return x

-- parsea una serie de expresiones separadas por espacios y luego aplica el constructor List.
parseList :: Parser LispVal
parseList =  List <$> sepBy parseExpr spaces

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

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

-- Funciones Auxiliares
-- parser que reconoce uno de los simbolos permitidos
-- en los identificadores de Scheme
symbol :: Parser Char
symbol = oneOf "!#$%&|*+-/:<=>?@^_~"

-- parser que reconoce espacios
spaces :: Parser ()
spaces = skipMany1 space

max_stars_repo_name: JuanBono/my-scheme

max_issues_repo_name: JuanBono/my-scheme

max_forks_repo_name: JuanBono/my-scheme


text:
{-|
Module      : Data.Meta.NMA
Description : Second stage NMA
Copyright   : (c) Thodoris Papakonstantinou, 2019
License     : GPL-3
Maintainer  : mail@tpapak.com
Stability   : experimental
Portability : POSIX

Given a study graph of direct effects perform network-meta analysis
-}

{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE NamedFieldPuns #-}

module Data.Meta.NMA
  ( NetworkEffects (..)
  , Spring (..)
  , springNMA
  )
where

import           Control.Applicative
import           Data.List
import Data.Tuple
import qualified Data.Map.Strict               as Map
import qualified Data.Set                      as Set
import           Data.Maybe
import           Data.Either
import qualified Data.IntMap.Strict as IM

import qualified Numeric.LinearAlgebra         as LA

import qualified Data.Graph.AdjacencyList      as G
import qualified Data.Graph.AdjacencyList.BFS      as BFS

import           Data.Meta.Effects
import Data.Meta.Matrices
import Data.Meta.Studies

data Effect a => NetworkEffects a =
  NetworkEffects { directs :: StudyGraph a
                 --, leagueTable :: LA.Matrix Double
                 --, lstar :: Map.Map G.Edge (Either Spring Spring)
                 --, lstar :: Map.Map G.Edge Double
                 , networkEstimates :: Map.Map TreatmentId (Map.Map TreatmentId Double)
                 , networkVariances :: Map.Map TreatmentId (Map.Map TreatmentId Double)
                 , debugmsg :: String
                 }
  deriving (Show, Eq)

-- | Matrices needed for solving the linear system
data Matrices = Matrices { kMatrix :: LA.Matrix Double
                         , llMatrix :: [Double]
                         , cMatrix :: LA.Matrix Double
                         , dMatrix :: LA.Matrix Double
                         , xMatrix :: LA.Matrix Double
                         , spnTr :: Set.Set G.Edge
                         , parents :: IM.IntMap G.Vertex
                         , vss :: [G.Vertex]
                         , directEdges :: Set.Set G.Edge
                         }

type NaturalLength = Double
type Hardness = Double
type SpringLength = Double

data Spring = Spring Hardness NaturalLength SpringLength
  deriving (Show, Eq)

reverseSpring :: Spring -> Spring
reverseSpring (Spring k ll l) = Spring k (-ll) (-l)

hardness :: Spring -> Hardness
hardness (Spring k _ _) = k

naturalLength :: Spring -> NaturalLength
naturalLength (Spring _ ll _) = ll

springLength :: Spring -> SpringLength
springLength (Spring _ _ l) = l

type TauSquare = Double

springNMA :: (ArmEffect a b, Gaussian a, Gaussian b) => [Study]
                             -> Maybe G.Vertex
                             -> Maybe TauSquare
                             -> (Arm -> Either String a)
                             -> Either String (NetworkEffects b)
springNMA studies mpinnedVertex mtau2 armEffect =
  let vs = foldl' (\acvs study ->
                    let sid = getStudyId study
                        arms = getStudyArms study
                     in case mtau2 of
                          Nothing -> Set.insert (Left sid) $ Set.union acvs
                            $ Set.fromList $ map (Right . tidOfArm) arms
                          Just tau2 ->
                            let tauvs = Set.fromList $ map (\a
                                      -> Left $ StudyId $ StringId ((show sid) ++ "_" ++ (show $ tidOfArm a))) arms
                                tidvs = Set.fromList $ map (Right . tidOfArm) arms
                             in Set.insert (Left sid) $ Set.union acvs $ Set.union tauvs tidvs
                  ) Set.empty studies  :: Set.Set (Either StudyId TreatmentId)
      vssts = IM.fromList $ zip [1..] $ Set.toList vs
      stsvs = Map.fromList $ map swap $ IM.toList vssts
      Right vsts = sequence $ IM.filter isRight vssts -- only treatment vertices
      tsvs = Map.fromList $ map swap $ IM.toList vsts
      esprings = foldl' (\eaces study ->
                    let sid = getStudyId study
                        arms = getStudyArms study
                        studysprings = sequence $  case armEffect $ head arms of
                          Left err -> [Left err]
                          Right aArm ->
                            foldl' (\acc a -> case armEffect a of
                                      Left err -> Left err : acc
                                      Right armef ->
                                        let tid = tidOfArm a
                                            tausid = StudyId $ StringId ((show sid) ++ "_" ++ show tid)
                                            vs = stsvs Map.! (Left sid)
                                            vt = stsvs Map.! (Right tid)
                                            vtau = stsvs Map.! (Left tausid)
                                            eff = relatedRelative armef
                                            nl = expectation eff
                                            k = 1 / (variance eff)
                                            spr = Spring k nl nl
                                         in  case mtau2 of
                                              Nothing ->
                                                let sttrEdge = G.Edge vs vt
                                                 in (Right (sttrEdge, spr)) : acc
                                              Just tau2 ->
                                                let sttauEdge = G.Edge vs vtau
                                                    sttauspr = Spring (2/tau2) 0 0 -- That's where tau2 is inserted
                                                    tautrEdge = G.Edge vtau vt
                                                 in [Right (sttauEdge, sttauspr)
                                                    , Right (tautrEdge, spr)] <> acc
                                     ) [] arms :: [Either String (G.Edge, Spring)]
                     in case eaces of
                          Left err -> Left err
                          Right aces ->
                             case studysprings of
                               Left err -> Left err
                               Right sprs -> Right (sprs ++ aces)
              ) (Right []) studies :: Either String [(G.Edge,Spring)]
    in case esprings of
       Left err -> Left err
       Right springs ->
         let springsMap = Map.fromList springs
             es = Map.keys springsMap
             numEdges = length es
             directedSG = G.graphFromEdges es
             undirectedSG = G.makeUndirected directedSG
             getSpring :: G.Edge -> Spring -- from undirected edge to spring
             getSpring e =
               case Map.lookup e springsMap of
                 Nothing -> reverseSpring $ springsMap Map.! (G.reverseEdge e)
                 Just sp -> sp
             pinnedV =
               case mpinnedVertex of
                 Just p -> p
                 Nothing -> 1
             -- | get the Matrices corresponding to defined pinned Vertex
             getMatrices :: G.Vertex -> Matrices
             getMatrices pinnedVertex =
             -- | Networks spanning tree
               let netbfs = BFS.bfs undirectedSG pinnedVertex
                   spnTr = Set.fromList $ BFS.spanningTree netbfs
                   parents = BFS.parent netbfs
                   vss = map (\(G.Edge u v)->v) (Set.toList spnTr)
                   llmap = Map.fromList $ map (\e ->
                              case Set.member e spnTr of
                                True -> (e, Right $ getSpring e)
                                False -> let re = G.reverseEdge e
                                          in case Set.member re spnTr of
                                               True -> (re, Right $ getSpring re)
                                               False -> (e, Left $ getSpring e)
                                        ) es
                   -- | L vector (E><1) with the natural lengths of springs which  respect
                   -- the spanning tree orientation
                   directEdges = Set.fromList $ Map.keys llmap
                   llMatrix = map (\de ->
                     let ms = llmap Map.! de
                      in case ms of
                          Left sp -> naturalLength sp
                          Right sp -> naturalLength sp)
                              $ Set.toList directEdges
                   -- weight matrix (E x E)
                   kMatrix = LA.diagl $ map (hardness . getSpring) $ Set.toList directEdges
                   -- cMatrix (V-1 x E)
                   cMatrix =
                     let vertexToRow :: G.Vertex -> [Double]
                         vertexToRow x = map (\(G.Edge u v) ->
                                         if x==u
                                            then 1.0
                                            else
                                             if x==v then (-1)
                                                     else 0
                                             )
                                       $ Set.toList directEdges
                      in LA.fromLists $ map vertexToRow vss
                   -- D matrix (V-1 x E)
                   dMatrix = cMatrix LA.<> kMatrix
                   -- X matrix (E x (V-1)) consistency matrix
                   -- Definition: X <> L* = L
                   xMatrix =
                     let edgeToRow :: G.Edge -> [Double]
                         edgeToRow (G.Edge u v) =
                           let r = replicate (length vss) 0
                               r' = IM.fromList $ zip vss r
                               followEnd :: G.Vertex
                                         -> IM.IntMap Double
                                         -> IM.IntMap Double
                               followEnd k r =
                                 if k == pinnedVertex
                                    then r
                                    else
                                      let r' = IM.adjust (+ 1) k r
                                          k' = (BFS.parent netbfs) IM.! k
                                       in followEnd k' r'
                               followStart k r =
                                 if k == pinnedVertex
                                    then r
                                    else
                                      let r' = IM.adjust (+ (-1)) k r
                                          k' = (BFS.parent netbfs) IM.! k
                                       in followStart k' r'
                               r'' = followEnd v $ followStart u r'
                            in map (\v -> r'' IM.! v) vss -- has to respect vss ordering
                      in LA.fromLists $ map edgeToRow $ Set.toList directEdges
                in Matrices { kMatrix
                            , llMatrix
                            , cMatrix
                            , dMatrix
                            , xMatrix
                            , parents
                            , spnTr
                            , vss
                            , directEdges
                            }
             -- | l* (V-1 >< 1) where V-1 are the end vertices of the spanning
             -- tree
             mtrcs = getMatrices pinnedV
             kMatrix'=  kMatrix mtrcs
             llMatrix'= llMatrix mtrcs
             cMatrix'=  cMatrix mtrcs
             dMatrix'=  dMatrix mtrcs
             xMatrix'=  xMatrix mtrcs
             parents'=  parents mtrcs
             spnTr'  =  spnTr mtrcs
             vss' = vss mtrcs
             directEdges' = directEdges mtrcs
             lstar = ((LA.inv (dMatrix' LA.<> xMatrix')) LA.<> dMatrix') LA.<> (numEdges LA.>< 1 $ llMatrix')
             lstarMap = Map.fromList $ zip (Set.toList spnTr') (LA.toList $ LA.flatten lstar)
             getLength :: G.Edge
                       -> G.Vertex -- pinned Vertex
                       -> Map.Map G.Edge Double -- lstarMap
                       -> IM.IntMap G.Vertex -- parents
                       -> NaturalLength
             -- | follows back the spanning tree adding parent edges' length
             getLength (G.Edge u v) pinnedV lstarMap parents =
               let followEnd :: G.Vertex
                       -> Double
                       -> Double
                   followEnd k l =
                       if k == pinnedV
                          then l
                          else
                            let k' = parents IM.! k
                                l' = lstarMap Map.! (G.Edge k' k)
                             in followEnd k' (l+l')
                   followStart k l =
                       if k == pinnedV
                          then l
                          else
                            let k' = parents IM.! k
                                l' = lstarMap Map.! (G.Edge k' k)
                             in followStart k' (l-l')
                in if u == v
                      then 0
                      else
                        followEnd v $ followStart u 0
             netests = foldl' (\acc utid ->
               let vmap = foldl'
                     (\mac vtid ->
                       let u = tsvs Map.! utid
                           v = tsvs Map.! vtid
                           l = getLength (G.Edge u v) pinnedV lstarMap parents'
                        in Map.insert vtid l mac
                     ) Map.empty $ Map.keys tsvs
                in Map.insert utid vmap acc
                   ) Map.empty $ Map.keys tsvs
             effectiveVariance :: G.Edge
                               -> G.Vertex -- pinned
                               -> [G.Vertex] -- vss
                               -> Set.Set G.Edge -- spanning tree
                               -> LA.Matrix Double -- D matrix
                               -> LA.Matrix Double -- X matrix
                               -> IM.IntMap G.Vertex -- parents
                               -> Hardness
             effectiveVariance (G.Edge u v) pinnedV vss spnTr dMatrix xMatrix parents =
               let r = IM.fromList $ zip vss $ replicate (length vss) 0
                   r' = case IM.lookup u r of
                               Nothing -> r
                               Just _ -> IM.insert u (-1.0) r
                   r'' = case IM.lookup v r' of
                               Nothing -> r'
                               Just _ -> IM.insert v (1.0) r'
                   forceVector = LA.tr' $ LA.fromLists [map (\v -> r'' IM.! v) vss] -- has to respect vss ordering
                   eflstar = (LA.inv (dMatrix LA.<> xMatrix)) LA.<> forceVector
                   eflstarMap = Map.fromList $
                     zip (Set.toList spnTr) (LA.toList $ LA.flatten eflstar)
                in abs $ getLength (G.Edge u v) pinnedV eflstarMap parents
             netvars = foldl' (\acc utid ->
               let vmap = foldl'
                     (\mac vtid ->
                       let u' = tsvs Map.! utid
                           v' = tsvs Map.! vtid
                           u = min u' v'
                           v = max u' v'
                           mtc = getMatrices u
                           l = if u == v
                                  then 0
                                  else
                                  --effectiveVariance (G.Edge u v) u (vss mtc) (spnTr mtc) (dMatrix mtc) (xMatrix mtc) (parents mtc)
                                  effectiveVariance (G.Edge u v) pinnedV (vss') (spnTr') (dMatrix') (xMatrix') (parents')
                        in Map.insert vtid l mac
                     ) Map.empty $ Map.keys tsvs
                in Map.insert utid vmap acc
                   ) Map.empty $ Map.keys tsvs
             sg = StudyGraph { studyGraph = G.graphFromEdges []
                             , directGraph = G.graphFromEdges $ Set.toList spnTr'
                             , vsts = vsts
                             , tsvs = tsvs
                             }
             nets = Map.empty
          in Right $ NetworkEffects { directs = sg
                                    --, leagueTable = dMatrix'
                                    --, lstar = lstar
                                    , networkEstimates = netests
                                    , networkVariances = netvars
                                    , debugmsg = unlines [ "spnTr " <> (show $ Set.toList spnTr')
                                                         , "vss " <> (show $ map (\(G.Edge u v)->v) (Set.toList spnTr'))
                                                         , "netvars" <> (show netvars)
                                                         , "dMatrix" <> (show dMatrix')
                                                         , "directEdges" <> (show directEdges')
                                                         ]
                                    }

max_stars_repo_name: tpapak/meta-analysis

max_issues_repo_name: tpapak/meta-analysis

max_forks_repo_name: tpapak/meta-analysis


text:
module Svg where

import Data.Maybe
import Control.Applicative
import Data.Complex

data ViewBox = ViewBox {
  viewBoxX :: String,
  viewBoxY :: String,
  viewBoxWidth :: String,
  viewBoxHeight :: String
  } deriving Eq

viewBoxToString :: ViewBox -> String
viewBoxToString vb = "viewBox=\"" ++ viewBoxX vb ++ " " ++ viewBoxY vb ++ " " ++
                     viewBoxWidth vb ++ " " ++ viewBoxHeight vb ++ "\""

data WidthHeight = WidthHeight {
  width :: String,
  height :: String
  } deriving Eq

switchOrientation :: WidthHeight -> WidthHeight
switchOrientation (WidthHeight w h) = WidthHeight h w

widthHeightToString :: WidthHeight -> String
widthHeightToString (WidthHeight w h) = "width=\"" ++ w ++ "\" height=\"" ++ h ++ "\""

a4PortraitWidthHeight :: WidthHeight
a4PortraitWidthHeight = WidthHeight "210mm" "297mm"

data Stroke = Stroke {
  strokeColorMaybe :: Maybe String,
  strokeWidthMaybe :: Maybe Double
  } deriving Eq

strokeToString :: Stroke -> String
strokeToString (Stroke cm wm) = maybeToAttr cm colorToString ++
                                maybeToAttr wm widthToString where
                                  colorToString c = "stroke=\"" ++ c ++ "\""
                                  widthToString w = "stroke-width=\"" ++ show w ++ "\""

strokeWidthMap :: (Double -> Double) -> Stroke -> Stroke
strokeWidthMap f (Stroke cm wm) = Stroke cm (f <$> wm)

maybeToAttr :: Maybe a -> (a -> String) -> String
maybeToAttr maybe toStr =
  fromMaybe "" (((' ' :) . toStr) <$> maybe)

maybePrefix :: Maybe String -> String
maybePrefix m = fromMaybe "" ((++ ":") <$> m)

prologue :: Maybe String -> String -> Maybe WidthHeight -> Maybe ViewBox -> String
prologue svgPrefixMaybe xlinkPrefix widthHeightMaybe viewBoxMaybe =
  "<" ++  maybePrefix svgPrefixMaybe  ++
  "svg xmlns" ++ (fromMaybe "" ((':' :) <$> svgPrefixMaybe)) ++
  "=\"http://www.w3.org/2000/svg\" xmlns:" ++
  xlinkPrefix ++ "=\"http://www.w3.org/1999/xlink\"" ++
  maybeToAttr widthHeightMaybe widthHeightToString ++
  maybeToAttr viewBoxMaybe viewBoxToString ++
  "><" ++ maybePrefix svgPrefixMaybe ++ "g transform=\"scale(1 -1)\">"

epilogue :: Maybe String -> String
epilogue svgPrefixMaybe = let pref = maybePrefix svgPrefixMaybe in
  "</" ++ pref ++ "g></" ++ pref ++ "svg>"

drawGraph :: Maybe String -> Double -> Double -> Stroke -> [(Double, Double)] -> String
drawGraph svgPrefixMaybe xScale yScale stroke pnts =
  let scaledPnts = map (\(x, y) -> (x * xScale, y * yScale)) pnts
      pairToString (x, y) = show x ++ " " ++ show y
  in
   "<" ++ maybePrefix svgPrefixMaybe ++ "path" ++
   strokeToString stroke ++
   " fill=\"none\" d=\"M " ++ pairToString (head scaledPnts) ++
   " L" ++ (foldr (\s1 s2 -> (' ' : s1) ++ s2) "" (map pairToString (tail scaledPnts))) ++
   "\" />"

drawComplex :: Show a => Complex a -> Maybe String -> Stroke -> String
drawComplex (a :+ b) svgPrefixMaybe stroke =
  "<" ++ maybePrefix svgPrefixMaybe ++ "line" ++ strokeToString stroke ++
  " x1=\"0\" y1=\"0\" x2=\"" ++ show a ++ "\" y2=\"" ++ show b ++ "\" />"

max_stars_repo_name: tilgalas/haskell-fourier

max_issues_repo_name: tilgalas/haskell-fourier

max_forks_repo_name: tilgalas/haskell-fourier


text:
{-# LANGUAGE FlexibleInstances, RankNTypes #-}
module Dna.Nucleotides (
    Nucleotide, Dna, DnaList, Rna, (+:+),
    isDnaPrefixOf, lengthDna,
    groupDna, tailDna, readsDna, nullDna,
    liftDna2, liftDna, takeDna, dnaToRna,
    dropDna, dnaComplement, FastaSequence,
    parseFastaFile, highestGC)
where

import           Text.Printf
import           Data.Ratio
import           Control.Monad.IO.Class
import           Data.List
import           Data.Functor
import           Data.STRef
import           Control.Monad
import           Control.Applicative
import           Numeric.LinearAlgebra
import           Numeric.LinearAlgebra.Data
import           Numeric.LinearAlgebra.Devel
import           Control.Monad.ST
import qualified Text.Parser.Combinators as P
import qualified Text.Parser.Char as P
import qualified Text.Parser.Token as P
import qualified Text.Trifecta.Parser as P

import qualified Data.ByteString.Lazy.Char8 as C
import qualified Data.Set as S

data Nucleotide = A | C | G | T | U
    deriving (Eq, Ord)

instance Show Nucleotide where
    show A = "A"
    show C = "C"
    show G = "G"
    show T = "T"
    show U = "U"

data FastaSequence a = Fasta {
    _fastaId :: String,
    _fastaInfo :: a}

class FastaParseInfo a where
    parse :: P.TokenParsing m => m a

instance Functor FastaSequence where
    fmap f (Fasta ident datum) = Fasta ident (f datum)

instance FastaParseInfo Dna where
    parse = DnaIntern . concat <$> P.sepEndBy (many parseDnaNucleotide) P.newline

instance FastaParseInfo Rna where
    parse = RnaIntern . concat <$> P.sepEndBy (many parseRnaNucleotide) P.newline

newtype Dna = DnaIntern [Nucleotide]
    deriving (Ord, Eq)
newtype Rna = RnaIntern [Nucleotide]
    deriving (Ord, Eq)
type DnaList = [Dna]

parseFastaFile :: (FastaParseInfo a, MonadIO m) => String -> m (Maybe [FastaSequence a])
parseFastaFile = P.parseFromFile parseFastaList

parseDnaNucleotide :: P.TokenParsing m => m Nucleotide
parseDnaNucleotide = P.choice [P.string "A" *> pure A, P.string "T" *> pure T, P.string "G" *> pure G, P.string "C" *> pure C]

parseRnaNucleotide :: P.TokenParsing m => m Nucleotide
parseRnaNucleotide = P.choice [P.string "A" *> pure A, P.string "U" *> pure U, P.string "G" *> pure G, P.string "C" *> pure C]

parseFastaList :: (FastaParseInfo a, P.TokenParsing m) => m [FastaSequence a]
parseFastaList = P.many parseFasta

parseFasta :: (FastaParseInfo a, P.TokenParsing m) => m (FastaSequence a)
parseFasta = Fasta
                <$> (P.string ">" *> P.manyTill P.anyChar P.newline)
                <*> parse

unDna (DnaIntern dna) = dna
unRna (DnaIntern rna) = rna

liftDna :: ([Nucleotide] -> [Nucleotide]) -> Dna -> Dna
liftDna f (DnaIntern dna) = DnaIntern (f dna)

liftDna2 :: ([Nucleotide] -> [Nucleotide] -> a) -> Dna -> Dna -> a
liftDna2 f (DnaIntern dna1) (DnaIntern dna2) =  f dna1 dna2

groupDna :: Dna -> DnaList
groupDna = map DnaIntern . group . unDna

dropDna :: Int -> Dna -> Dna
dropDna n = liftDna $ drop n

takeDna :: Int -> Dna -> Dna
takeDna n = liftDna $ take n

tailDna :: Dna -> Dna
tailDna = liftDna tail

lengthDna :: Dna -> Int
lengthDna (DnaIntern dna) = length dna

isDnaPrefixOf :: Dna -> Dna -> Bool
isDnaPrefixOf  = liftDna2 isPrefixOf

(+:+) :: Dna -> Dna -> Dna
l +:+ r = DnaIntern $ liftDna2 (++) l r

nullDna :: Dna
nullDna = DnaIntern []

instance Show Dna where
    show (DnaIntern dna) = foldr ((++) . show) "" dna

instance Show Rna where
    show (RnaIntern rna) = foldr ((++) . show) "" rna

instance (Show a) => Show (FastaSequence a) where
    show (Fasta ident datum) = ">" ++ ident ++ "\n" ++ show datum ++ "\n"

instance {-# OVERLAPPING #-} Show (FastaSequence Rational) where
    show (Fasta ident datum) = ">" ++ ident ++ "\n" ++ show (fromRational datum :: Double)  ++ "\n"

instance {-# OVERLAPPING #-} Show DnaList where
    show = intercalate "\n" . map show

readsNucleotide          :: Int -> ReadS Nucleotide
readsNucleotide _ ('A':xs) = [(A,xs)]
readsNucleotide _ ('C':xs) = [(C,xs)]
readsNucleotide _ ('G':xs) = [(G,xs)]
readsNucleotide _ ('T':xs) = [(T,xs)]
readsNucleotide _ ('U':xs) = [(U,xs)]

instance Read Nucleotide where
    readsPrec = readsNucleotide

highestGC :: [FastaSequence Dna] -> FastaSequence Rational
highestGC = last . sortOn  _fastaInfo . fmap  ((`frequency` [C,G]) <$>)

readsDna        :: Int -> ReadS Dna
readsDna  _ []  = [(DnaIntern [],[])]
readsDna  d s@(x:xs)
    | x `elem` ['A','C','G','T'] = [(DnaIntern $ read [x]:dna, rest) | (DnaIntern dna,rest) <- readsDna (d+1) xs ]
    | otherwise = [(DnaIntern [],s)]

readsRna        :: Int -> ReadS Rna
readsRna  _ []  = [(RnaIntern [],[])]
readsRna  d s@(x:xs)
    | x `elem` ['A','C','G','U'] = [(RnaIntern $ read [x]:rna, rest) | (RnaIntern rna,rest) <- readsRna (d+1) xs ]
    | otherwise = [(RnaIntern [],s)]

readsDnaList    :: Int -> ReadS DnaList
readsDnaList _ [] = [([],[])]
readsDnaList d ('\n':xs) = readsDnaList (d+1) xs
readsDnaList d s@(x:xs)
    | x `notElem` ['A','C','G','T'] = [([],s)]
    | otherwise = [(DnaIntern dna:restDna,rest) |
            (DnaIntern dna,restRead) <- readsDna d s,
            (restDna,rest) <- readsDnaList (d+length dna) restRead]

mapComplementDan A = T
mapComplementDan T = A
mapComplementDan C = G
mapComplementDan G = C

dnaComplement :: Dna -> Dna
dnaComplement = liftDna $ reverse . map mapComplementDan

instance Read Dna where
    readsPrec = readsDna

instance Read Rna where
    readsPrec = readsRna

instance {-# OVERLAPPING #-} Read DnaList where
    readsPrec = readsDnaList

isConnected :: Dna -> Dna -> Bool
isConnected (DnaIntern first) (DnaIntern second) = and $ zipWith (==) (tail first) second

dnaToRna :: Dna -> Rna
dnaToRna = RnaIntern . map mapNucleotidDnaRna . unDna
    where mapNucleotidDnaRna T = U
          mapNucleotidDnaRna a = a

frequency :: Dna -> [Nucleotide] -> Rational
frequency (DnaIntern d) pop = fromIntegral (length (filter (`elem` pop) d)) %  fromIntegral (length d)

max_stars_repo_name: mathk/NautilusDna

max_issues_repo_name: mathk/NautilusDna

max_forks_repo_name: mathk/NautilusDna


text:
--
-- FlattenLayer: flatten/unflatten layer
--

module CNN.FlattenLayer (
  flatten
, unflatten
) where

import qualified Numeric.LinearAlgebra as A

import CNN.Image
import CNN.LayerType

{- |
flatten: flatten image

  IN: Image (=[Matrix R])

>>> let im = [A.fromLists [[2.0,3.0],[4.0,5.0]], A.fromLists [[1.0,2.0],[4.0,3.0]], A.fromLists [[2.0,1.0],[6.0,5.0]]]
>>> flatten im
[(1><12)
 [ 2.0, 3.0, 4.0, 5.0, 1.0, 2.0, 4.0, 3.0, 2.0, 1.0, 6.0, 5.0 ]]

-}

flatten :: Image -> Image
flatten im = [A.fromLists [concat $ map (A.toList . A.flatten) im]]

{- |
unflatten: unflatten image

  IN: [Matrix R]

>>> let ds = [A.fromLists [[2.0,3.0,4.0,5.0,1.0,2.0,4.0,3.0,2.0,1.0,6.0,5.0]]]
>>> unflatten 2 2 ds
[(2><2)
 [ 2.0, 3.0
 , 4.0, 5.0 ],(2><2)
 [ 1.0, 2.0
 , 4.0, 3.0 ],(2><2)
 [ 2.0, 1.0
 , 6.0, 5.0 ]]

-}

unflatten :: Int -> Int -> Image -> Image
unflatten x y [m] = map (A.reshape x . A.fromList) $ split (y*x) v
  where
    v = A.toList $ A.flatten m
    split :: Int -> [a] -> [[a]]
    split _ [] = []
    split n ds = d : split n ds'
      where
        (d, ds') = splitAt n ds

max_stars_repo_name: eijian/deeplearning

max_issues_repo_name: eijian/deeplearning

max_forks_repo_name: eijian/deeplearning


text:
{-# LANGUAGE CPP                #-}
{-# LANGUAGE DeriveGeneric      #-}
{-# LANGUAGE FlexibleInstances  #-}
{-# LANGUAGE StandaloneDeriving #-}

-- | Flat instances for the base library
module Flat.Instances.Base () where

import           Control.Monad         (liftM2)
import           Data.Bool
import           Data.Char
import           Data.Complex          (Complex (..))
import           Data.Fixed
-- #if MIN_VERSION_base(4,9,0)
import qualified Data.List.NonEmpty    as B
-- #endif

#if ! MIN_VERSION_base(4,8,0)
import           Control.Applicative
import           Data.Monoid           (mempty)
#endif

#if MIN_VERSION_base(4,9,0)
import qualified Data.Semigroup        as Semigroup
#endif

import qualified Data.Monoid           as Monoid
import           Data.Ratio
import           Flat.Instances.Util
import           Prelude               hiding (mempty)

#if !MIN_VERSION_base(4,11,0)
import           Data.Monoid           ((<>))
#endif

#if MIN_VERSION_base(4,8,0)
import           Data.Functor.Identity (Identity (..))
#endif

#if !MIN_VERSION_base(4,9,0)
deriving instance Generic (Complex a)
#endif

{- ORMOLU_DISABLE -}
-- $setup
-- >>> :set -XNegativeLiterals
-- >>> import Flat.Instances.Test
-- >>> import Data.Fixed
-- >>> import Data.Int
-- >>> import Data.Complex(Complex(..))
-- >>> import Numeric.Natural
-- >>> import Data.Word
-- >>> import Data.Ratio
-- >>> import Flat.Run
-- >>> import Data.Monoid
-- >>> import qualified Data.List.NonEmpty as B
-- >>> let test = tstBits
-- >>> let y = 33
{- ORMOLU_ENABLE -}

-- >>> y

#if MIN_VERSION_base(4,8,0)
-- | @since 0.4.4
instance Flat a => Flat (Identity a) where
    encode (Identity a) = encode a
    size (Identity a) = size a
    decode = Identity <$> decode
#endif

-- | @since 0.4.4
instance Flat a => Flat (Monoid.Dual a) where
    encode (Monoid.Dual a) = encode a
    size (Monoid.Dual a) = size a
    decode = Monoid.Dual <$> decode

-- | @since 0.4.4
instance Flat Monoid.All where
    encode (Monoid.All a) = encode a
    size (Monoid.All a) = size a
    decode = Monoid.All <$> decode

-- | @since 0.4.4
instance Flat Monoid.Any where
    encode (Monoid.Any a) = encode a
    size (Monoid.Any a) = size a
    decode = Monoid.Any <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Monoid.Sum a) where
    encode (Monoid.Sum a) = encode a
    size (Monoid.Sum a) = size a
    decode = Monoid.Sum <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Monoid.Product a) where
    encode (Monoid.Product a) = encode a
    size (Monoid.Product a) = size a
    decode = Monoid.Product <$> decode

#if MIN_VERSION_base(4,9,0)
{- |

>>> let w = Just (11::Word8); a = Alt w <> Alt (Just 24) in tst a == tst w
True

>>> let w = Just (11::Word8); a = Alt Nothing <> Alt w in tst a == tst w
True

@since 0.4.4
-}
instance Flat (f a) => Flat (Monoid.Alt f a) where
    encode (Monoid.Alt a) = encode a
    size (Monoid.Alt a) = size a
    decode = Monoid.Alt <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Semigroup.Min a) where
    encode (Semigroup.Min a) = encode a
    size (Semigroup.Min a) = size a
    decode = Semigroup.Min <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Semigroup.Max a) where
    encode (Semigroup.Max a) = encode a
    size (Semigroup.Max a) = size a
    decode = Semigroup.Max <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Semigroup.First a) where
    encode (Semigroup.First a) = encode a
    size (Semigroup.First a) = size a
    decode = Semigroup.First <$> decode

-- | @since 0.4.4
instance Flat a => Flat (Semigroup.Last a) where
    encode (Semigroup.Last a) = encode a
    size (Semigroup.Last a) = size a
    decode = Semigroup.Last <$> decode
#endif

{- |
`()`, as all data types with a single constructor, has a zero-length encoding.

>>> test ()
(True,0,"")
-}
instance Flat () where
    encode _ = mempty

    size _ = id

    decode = pure ()

{- |
One bit is plenty for a Bool.

>>> test False
(True,1,"0")

>>> test True
(True,1,"1")
-}
instance Flat Bool where
    encode = eBool

    size = sBool

    decode = dBool

{- |
Char's are mapped to Word32 and then encoded.

For ascii characters, the encoding is standard ascii.

>>> test 'a'
(True,8,"01100001")

For unicode characters, the encoding is non standard.

>>> test 'È'
(True,16,"11001000 00000001")

>>> test '不'
(True,24,"10001101 10011100 00000001")

#ifndef ETA
>>> test "\x1F600"
(True,26,"11000000 01110110 00000011 10")
#endif
-}
instance Flat Char where
    size = sChar

    encode = eChar

    decode = dChar

{- |
>>> test (Nothing::Maybe Bool)
(True,1,"0")

>>> test (Just False::Maybe Bool)
(True,2,"10")
-}
instance Flat a => Flat (Maybe a)

{- |
>>> test (Left False::Either Bool ())
(True,2,"00")

>>> test (Right ()::Either Bool ())
(True,1,"1")
-}
instance (Flat a, Flat b) => Flat (Either a b)

{- |
>>> test (MkFixed 123 :: Fixed E0)
(True,16,"11110110 00000001")

>>> test (MkFixed 123 :: Fixed E0) == test (MkFixed 123 :: Fixed E2)
True
-}
instance Flat (Fixed a) where
    encode (MkFixed n) = encode n

    size (MkFixed n) = size n

    decode = MkFixed <$> decode

{- |
Word8 always take 8 bits.

>>> test (0::Word8)
(True,8,"00000000")

>>> test (255::Word8)
(True,8,"11111111")
-}
instance Flat Word8 where
    encode = eWord8

    decode = dWord8

    size = sWord8

{- |
Natural, Word, Word16, Word32 and Word64 are encoded as a non empty list of 7 bits chunks (least significant chunk first and most significant bit first in every chunk).

Words are always encoded in a whole number of bytes, as every chunk is 8 bits long (1 bit for the List constructor, plus 7 bits for the value).

The actual definition is:

@
Word64 ≡   Word64 Word

Word32 ≡   Word32 Word

Word16 ≡   Word16 Word

Word ≡   Word (LeastSignificantFirst (NonEmptyList (MostSignificantFirst Word7)))

LeastSignificantFirst a ≡   LeastSignificantFirst a

NonEmptyList a ≡   Elem a
                 | Cons a (NonEmptyList a)

MostSignificantFirst a ≡   MostSignificantFirst a

Word7 ≡   V0
        | V1
        | V2
        ...
        | V127
@

Values between as 0 and 127 fit in a single byte.

127 (0b1111111) is represented as Elem V127 and encoded as: Elem=0 127=1111111

>>> test (127::Word)
(True,8,"01111111")

254 (0b11111110) is represented as Cons V126 (Elem V1) (254=128+126) and encoded as: Cons=1 V126=1111110 (Elem=0 V1=0000001):

>>> test (254::Word)
(True,16,"11111110 00000001")

Another example, 32768 (Ob1000000000000000 = 0000010 0000000 0000000):

>>> test (32768::Word32)
(True,24,"10000000 10000000 00000010")

As this is a variable length encoding, values are encoded in the same way, whatever their type:

>>> all (test (3::Word) ==) [test (3::Word16),test (3::Word32),test (3::Word64)]
True
-}
instance Flat Word where
    size = sWord

    encode = eWord

    decode = dWord

{- |
Naturals are encoded just as the fixed size Words.

>>> test (0::Natural)
(True,8,"00000000")

>>> test (2^120::Natural)
(True,144,"10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 00000010")
-}
instance Flat Natural where
    size = sNatural

    encode = eNatural

    decode = dNatural

instance Flat Word16 where
    encode = eWord16

    decode = dWord16

    size = sWord16

instance Flat Word32 where
    encode = eWord32

    decode = dWord32

    size = sWord32

instance Flat Word64 where
    encode = eWord64

    decode = dWord64

    size = sWord64

{- |
Integer, Int, Int16, Int32 and Int64 are defined as the <https://developers.google.com/protocol-buffers/docs/encoding#signed-integers ZigZag> encoded version of the equivalent unsigned Word:

@
Int   ≡  Int   (ZigZag Word)

Int64 ≡  Int64 (ZigZag Word64)

Int32 ≡  Int32 (ZigZag Word32)

Int16 ≡  Int16 (ZigZag Word16)

Int8  ≡  Int8  (ZigZag Word8)

ZigZag a ≡ ZigZag a
@

ZigZag encoding alternates between positive and negative numbers, so that numbers whose absolute value is small can be encoded efficiently:

>>> test (0::Int)
(True,8,"00000000")

>>> test (-1::Int)
(True,8,"00000001")

>>> test (1::Int)
(True,8,"00000010")

>>> test (-2::Int)
(True,8,"00000011")

>>> test (2::Int)
(True,8,"00000100")
-}
instance Flat Int where
    size = sInt

    encode = eInt

    decode = dInt

{- |
Integers are encoded just as the fixed size Ints.

>>> test (0::Integer)
(True,8,"00000000")

>>> test (-1::Integer)
(True,8,"00000001")

>>> test (1::Integer)
(True,8,"00000010")

>>> test (-(2^4)::Integer)
(True,8,"00011111")

>>> test (2^4::Integer)
(True,8,"00100000")

>>> test (-(2^120)::Integer)
(True,144,"11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000011")

>>> test (2^120::Integer)
(True,144,"10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 10000000 00000100")
-}
instance Flat Integer where
    size = sInteger

    encode = eInteger

    decode = dInteger

{- |
>>> test (0::Int8)
(True,8,"00000000")

>>> test (127::Int8)
(True,8,"11111110")

>>> test (-128::Int8)
(True,8,"11111111")
-}
instance Flat Int8 where
    encode = eInt8

    decode = dInt8

    size = sInt8

{- |
>>> test (0::Int16)
(True,8,"00000000")

>>> test (1::Int16)
(True,8,"00000010")

>>> test (-1::Int16)
(True,8,"00000001")

>>> test (minBound::Int16)
(True,24,"11111111 11111111 00000011")

equivalent to 0b1111111111111111

>>> test (maxBound::Int16)
(True,24,"11111110 11111111 00000011")

equivalent to 0b1111111111111110
-}
instance Flat Int16 where
    size = sInt16

    encode = eInt16

    decode = dInt16

{- |
>>> test (0::Int32)
(True,8,"00000000")

>>> test (minBound::Int32)
(True,40,"11111111 11111111 11111111 11111111 00001111")

>>> test (maxBound::Int32)
(True,40,"11111110 11111111 11111111 11111111 00001111")
-}
instance Flat Int32 where
    size = sInt32

    encode = eInt32

    decode = dInt32

{- |
>>> test (0::Int64)
(True,8,"00000000")

>>> test (minBound::Int64)
(True,80,"11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000001")

>>> test (maxBound::Int64)
(True,80,"11111110 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111 00000001")
-}
instance Flat Int64 where
    size = sInt64

    encode = eInt64

    decode = dInt64

{- |
Floats are encoded as standard IEEE binary32 values:

@
IEEE_754_binary32 ≡ IEEE_754_binary32 {sign :: Sign,
                                        exponent :: MostSignificantFirst Bits8,
                                        fraction :: MostSignificantFirst Bits23}
@

>>> test (0::Float)
(True,32,"00000000 00000000 00000000 00000000")

>>> test (1.4012984643E-45::Float)
(True,32,"00000000 00000000 00000000 00000001")

>>> test (1.1754942107E-38::Float)
(True,32,"00000000 01111111 11111111 11111111")
-}
instance Flat Float where
    size = sFloat

    encode = eFloat

    decode = dFloat

{- |
Doubles are encoded as standard IEEE binary64 values:

@
IEEE_754_binary64 ≡ IEEE_754_binary64 {sign :: Sign,
                                        exponent :: MostSignificantFirst Bits11,
                                        fraction :: MostSignificantFirst Bits52}
@
-}
instance Flat Double where
    size = sDouble

    encode = eDouble

    decode = dDouble

{- |
>>> test (4 :+ 2 :: Complex Word8)
(True,16,"00000100 00000010")
-}
instance Flat a => Flat (Complex a)

{- |
Ratios are encoded as tuples of (numerator,denominator)

>>> test (3%4::Ratio Word8)
(True,16,"00000011 00000100")
-}
instance (Integral a, Flat a) => Flat (Ratio a) where
    size a = size (numerator a, denominator a)

    encode a = encode (numerator a, denominator a)

    -- decode = uncurry (%) <$> decode
    decode = liftM2 (%) decode decode

{- |
>>> test ([]::[Bool])
(True,1,"0")

>>> test [False,False]
(True,5,"10100")

This instance and other similar ones are declared as @OVERLAPPABLE@, because for better encoding/decoding
performance it can be useful to declare instances of concrete types, such as @[Char]@ (not provided out of the box).
-}
instance {-# OVERLAPPABLE #-} Flat a => Flat [a]

{-
>>> import Weigh
>>> flat [1..10::Int]
-}

-- Generic list instance (stack overflows with ETA, see https://github.com/typelead/eta/issues/901)
-- where
--size [] n = n+1
--size (h:t) n = trampoline size t (trampoline size h (n+1))
-- size = sizeListWith size -- foldl' (\n e -> ) n
-- encode = error "BAD"
-- encode = trampoline . encodeListWith encode
-- decode = decodeListWith decode
-- sizeListWith siz l n = foldl' (\n e -> 1 + n + siz e 0) n l
-- #ifdef ETA_VERSION
-- import Data.Function(trampoline)
-- import GHC.IO(trampolineIO)
-- #else
-- trampoline = id
-- trampolineIO = id
-- #endif

-- #if MIN_VERSION_base(4,9,0)

{- |
>>> test (B.fromList [True])
(True,2,"10")

>>> test (B.fromList [False,False])
(True,4,"0100")
-}
instance {-# OVERLAPPABLE #-} Flat a => Flat (B.NonEmpty a)

-- #endif

{- |
Tuples are supported up to 7 elements.

>>> test (False,())
(True,1,"0")

>>> test ((),())
(True,0,"")

"7 elements tuples ought to be enough for anybody" (Bill Gates - apocryphal)

>>> test (False,True,True,True,False,True,True)
(True,7,"0111011")

tst (1::Int,"2","3","4","5","6","7","8")
...error
-}

-- Not sure if these should be OVERLAPPABLE
instance {-# OVERLAPPABLE #-} (Flat a, Flat b) => Flat (a, b)

instance {-# OVERLAPPABLE #-} (Flat a, Flat b, Flat c) => Flat (a, b, c)

instance
    {-# OVERLAPPABLE #-}
    (Flat a, Flat b, Flat c, Flat d) =>
    Flat (a, b, c, d)

instance
    {-# OVERLAPPABLE #-}
    (Flat a, Flat b, Flat c, Flat d, Flat e) =>
    Flat (a, b, c, d, e)

instance
    {-# OVERLAPPABLE #-}
    (Flat a, Flat b, Flat c, Flat d, Flat e, Flat f) =>
    Flat (a, b, c, d, e, f)

instance
    {-# OVERLAPPABLE #-}
    ( Flat a
    , Flat b
    , Flat c
    , Flat d
    , Flat e
    , Flat f
    , Flat g
    ) =>
    Flat (a, b, c, d, e, f, g)

max_stars_repo_name: Bodigrim/flat

max_issues_repo_name: Bodigrim/flat

max_forks_repo_name: Bodigrim/flat


text:
{-# LANGUAGE CPP                       #-}
{-# LANGUAGE BangPatterns              #-}
{-# LANGUAGE DataKinds                 #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE GADTs                     #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TypeFamilies              #-}

import           Control.Applicative
import           Control.Monad.Random
import           Options.Applicative

import           Data.Function                (on)
import           Data.List

import           Graphics.Image               hiding (map, on)

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

import           Grenade
import           Grenade.Utils.ImageNet

data ResNetOptions = ResNetOptions FilePath         -- onnx file
                                   FilePath         -- input image

pResnet :: Parser ResNetOptions
pResnet = ResNetOptions <$> argument str (metavar "onnx") <*> argument str (metavar "image")

-- loadResNetImage :: FilePath -> IO (Maybe (S ('D3 224 224 3)))
loadResNetImage :: FilePath -> IO (Maybe (S ('D3 224 224 3)))
loadResNetImage path = do
  img <- readImageRGB VU path
  displayImage img
  return $ do
    guard $ dims img == (224, 224)
    let [img_red, img_green, img_blue] = toImagesX img
        [reds, greens, blues]          = map (map (\(PixelX y) -> doubleToRealNum y) . concat . toLists) [img_red, img_green, img_blue]

        redM   = H.dmmap (\a -> (a - 0.485) / 0.229) . H.tr $ H.fromList reds   :: H.L 224 224
        greenM = H.dmmap (\a -> (a - 0.456) / 0.224) . H.tr $ H.fromList greens :: H.L 224 224
        blueM  = H.dmmap (\a -> (a - 0.406) / 0.225) . H.tr $ H.fromList blues  :: H.L 224 224
        mat    = redM H.=== greenM H.=== blueM

    return (S3D mat)

main :: IO ()
main = do
    ResNetOptions netPath imgPath <- execParser (info (pResnet <**> helper) idm)
    res <- loadResNet netPath

    inputM <- loadResNetImage imgPath

    case (res, inputM) of
      (Right net, Just input)  -> do
        let S1D y = runNet net input
            tops  = getTop 5 $ LA.toList $ H.extract y
        mapM_ (\(i :: Int, x) -> print $ show i ++ ": " ++ (show . getLabel) x) $ zip [1..] tops
      (Left err, _) -> print err
      _             -> error "Could not load image"
  where
    getTop :: Ord a => Int -> [a] -> [Int]
    getTop n xs = map fst $ take n $ sortBy (flip compare `on` snd) $ zip [0..] xs