text:
{-# LANGUAGE FlexibleInstances, MultiParamTypeClasses #-}

module Physics.Quantum.Examples.HarmonicOscillator where

import Physics.Quantum.Algebra

import Numeric.Additive.Class          (Additive((+)), Abelian)
import Numeric.Algebra.Class           ( Multiplicative((*), pow1p)
                                       , Monoidal(..), Semiring, LeftModule(..), RightModule(..)
                                       )

import Numeric.Algebra.Division        (Division(..))
import Numeric.Domain.Euclidean        (Euclidean(..))
import Numeric.Domain.PID              (PID(..))
import Numeric.Domain.UFD              (UFD)
import Numeric.Domain.GCD              (GCDDomain(..))
import Numeric.Domain.Integral         (IntegralDomain(..))
import Numeric.Algebra.Unital          (Unital(..))
import Numeric.Algebra.Unital.UnitNormalForm (UnitNormalForm(..))
import Numeric.Decidable.Zero          (DecidableZero(..))
import Numeric.Semiring.ZeroProduct    (ZeroProductSemiring)
import Numeric.Algebra.Commutative     (Commutative)
import Numeric.Decidable.Units         (DecidableUnits(..))
import Numeric.Ring.Class              (Ring(..))
import Numeric.Decidable.Associates    (DecidableAssociates(..))
import Numeric.Rig.Class               (Rig(..))
import Numeric.Additive.Group          (Group(..))

import GHC.Natural                     (Natural)
import GHC.Float                       (double2Float)

import Data.Array.Repa.Eval            (Elt(touch))

import Data.Complex                    (Complex(..), magnitude, realPart)

import Prelude hiding ((+),(-),(/),(*), (^))
import qualified Prelude as P

import qualified Data.Map as M
import qualified Data.Array.Repa as R
import Data.Array.Repa                              (Z(..), (:.)(..))
import Data.Functor.Identity (runIdentity)

import qualified Data.Eigen.SparseMatrix as ES
import qualified Data.Eigen.Matrix as ESD
import Data.Eigen.Matrix                            (CComplex)
import Foreign.C.Types                              (CDouble)

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

import Graphics.Gloss
import Graphics.Gloss.Data.ViewPort                 (ViewPort)
import Graphics.Gloss.Raster.Field                  (rgb', makePicture)

import Debug.Trace           (trace)
import Control.Concurrent    (threadDelay)

--------------------
-- Add dimensions to quantities
-- Since we are writing operators as endomorphisms of vector spaces
-- if we have to give dimensiones to operators/eigenvalues,
-- we may need to declare the algebra of dimensions as instance of field too
-- Or we may need to relax the condition of operator to homomorphims of vector spaces
-------------------

------------------------------------------
-- Configure
------------------------------------------

data OscillatorInABoxConfig = OscillatorInABoxConfig { plankConstantConfig :: Double
                                                     , particleMassConfig :: Double
                                                     , oscillatorFrequencyConfig :: Double
                                                     , halfWidthConfig :: Int
                                                     , halfHeightConfig :: Int
                                                     , initialAmplitudeConfig :: Position -> Complex Double
                                                     , fpsConfig :: Int
                                                     -- the simulation will play this times faster than real time
                                                     -- mainly to keep time steps small while keeping fps low
                                                     -- the time step in each simulation will be this divided by fps
                                                     , timeFactorConfig :: TimeStep
                                                     }

type TimeStep = Double

data Position = Position { posX :: Int , posY :: Int }
    deriving (Eq, Ord, Show)

type Variance = Double

gaussianWave :: Int -> Int -> Variance -> (Position -> Complex Double)
gaussianWave xCenter yCenter var = \(Position x y) ->
    ( exp ( (
              - ( ((fromIntegral (x - xCenter))P.^2) P.+ ((fromIntegral (y - yCenter) )P.^2) )
            ) P./ (2 P.* (var P.^ 2))
          )
    ) :+ 0

smallBoxConfig :: OscillatorInABoxConfig
smallBoxConfig = OscillatorInABoxConfig { plankConstantConfig = 6.62607004e-16 -- scaling length dimension by 10^9 so that 1 corresponds to 1 nm
                                        , particleMassConfig = 10e-31 -- roughly mass of an electron
                                        , oscillatorFrequencyConfig = 10000
                                        , halfWidthConfig = 10
                                        , halfHeightConfig = 10
                                          -- tune variance of wave based on mass and frequency (and Plank's constant) to get a nice simulation
                                        , initialAmplitudeConfig = gaussianWave 0 0 3.8845e-6
                                        , fpsConfig = 30
                                        , timeFactorConfig = 1
                                        }

bigBoxConfig :: OscillatorInABoxConfig
bigBoxConfig = OscillatorInABoxConfig { plankConstantConfig = 6.62607004e-16 -- scaling length dimension by 10^9 so that 1 corresponds to 1 nm
                                      , particleMassConfig = 10e-31
                                      , oscillatorFrequencyConfig = 10000
                                      , halfWidthConfig = 100
                                      , halfHeightConfig = 100
                                        -- tune variance of wave based on mass and frequency (and Plank's constant) to get a nice simulation
                                      , initialAmplitudeConfig = gaussianWave 0 0 3.8845e-6
                                      , fpsConfig = 30
                                      , timeFactorConfig = 1
                                      }

---------------------
-- Helper functions to generate lattice and translate between two dimensional position and one dimensional array
-- may not need to modify these
--------------------

positionRangeConfig :: OscillatorInABoxConfig -> [Position]
positionRangeConfig config = [Position x y | x <- [-halfWidth .. halfWidth], y <- [-halfHeight .. halfHeight]]
    where halfWidth = halfWidthConfig config
          halfHeight = halfHeightConfig config

indexRangeConfig :: OscillatorInABoxConfig -> Int
indexRangeConfig config = (2 * halfWidth + 1) * (2 * halfHeight + 1)
    where halfWidth = halfWidthConfig config
          halfHeight = halfHeightConfig config

positionToIndexConfig :: OscillatorInABoxConfig -> Position -> Int
positionToIndexConfig config pos = ((posX pos + halfWidth) * (2 * halfHeight + 1)) + (posY pos + halfHeight)
    where halfWidth = halfWidthConfig config
          halfHeight = halfHeightConfig config

indexToPositionConfig :: OscillatorInABoxConfig -> Int -> Position
indexToPositionConfig config index = Position ( (floor ((fromIntegral index) P./ (fromIntegral $! 2 * halfHeight + 1))) - halfWidth )
                                              ( (index `mod` (2 * halfHeight + 1)) - halfHeight )
    where halfWidth = halfWidthConfig config
          halfHeight = halfHeightConfig config

------------------------------------
-- simulate oscillation using Eigen library's sparse matrices
------------------------------------

sparseMatModelSimulator :: OscillatorInABoxConfig
                        -> ( TimeStep -> SparseMatOperator (Complex Double) (CComplex CDouble) Position
                           , SparseMatKet (Complex Double) (CComplex CDouble) Position
                           , SparseMatKet (Complex Double) (CComplex CDouble) Position -> Picture
                           , OscillatorInABoxConfig
                           )
sparseMatModelSimulator config = (oscillationOperator, initialModel, drawing, config)
    where
        -- shorthands for config params
        positionRange = positionRangeConfig config
        halfWidth = halfWidthConfig config
        halfHeight = halfHeightConfig config
        indexRange = indexRangeConfig config

        positionToIndex = positionToIndexConfig config
        indexToPosition = indexToPositionConfig config

        plankConstant = plankConstantConfig config
        particleMass = particleMassConfig config
        oscillatorFrequency = oscillatorFrequencyConfig config

        initialAmplitude = initialAmplitudeConfig config

        -- position and momentum operators

        posXOp :: SparseMatOperator (Complex Double) (CComplex CDouble) Position
        posXOp = SparseMatOperator $! ES.fromList indexRange indexRange $! map (\pos@(Position x _) -> (positionToIndex pos, positionToIndex pos, fromIntegral x)) positionRange

        posYOp :: SparseMatOperator (Complex Double) (CComplex CDouble) Position
        posYOp = SparseMatOperator $! ES.fromList indexRange indexRange $! map (\pos@(Position _ y) -> (positionToIndex pos, positionToIndex pos, fromIntegral y)) positionRange

        momXOp :: SparseMatOperator (Complex Double) (CComplex CDouble) Position
        -- taking dx as 2
        -- <x'| P_x |a>  =  -ih ( d/dx(a)|_x' ) = - ( a(x'+1) - a(x'-1) ) * ih/2   = - (<x'+1|a> - <x'-1|a>) * ih / 2
        momXOp = SparseMatOperator $! ES.fromList indexRange indexRange $! ( (map (\pos@(Position x y) -> ( positionToIndex pos
                                                                                , positionToIndex (Position (x P.+ 1) y)
                                                                                , 0 :+ (- plankConstant/2)
                                                                                )
                                                                                )
                                                                              -- removing max x from position range to keep operator within lattice
                                                                              -- that would give 0 contribution anyway
                                                                              [Position i j | i <- [(-halfWidth) .. (halfWidth-1)], j <- [-(halfHeight) .. halfHeight]]
                                                                             )
                                                                             ++
                                                                             (map (\pos@(Position x y) -> ( positionToIndex pos
                                                                                , positionToIndex (Position (x P.- 1) y)
                                                                                , 0 :+ (  plankConstant/2)
                                                                                )
                                                                                )
                                                                              -- removing min x from position range to keep operator within lattice
                                                                              -- that would give 0 contribution anyway
                                                                              [Position i j | i <- [-(halfWidth-1) .. (halfWidth)], j <- [-(halfHeight) .. halfHeight]]
                                                                             )
                                                                           )

        momYOp :: SparseMatOperator (Complex Double) (CComplex CDouble) Position
        -- taking dy as 2
        -- <y'| P_y |a>  =  -ih ( d/dy(a)|_y' ) = - ( a(y'+1) - a(y'-1) ) * ih/2   = - (<y'+1|a> - <y'-1|a>) * ih / 2
        momYOp = SparseMatOperator $! ES.fromList indexRange indexRange $! ( (map (\pos@(Position x y) -> ( positionToIndex pos
                                                                                , positionToIndex (Position x (y P.+ 1))
                                                                                , 0 :+ (- plankConstant/2)
                                                                                )
                                                                                )
                                                                              -- removing max y from position range to keep operator within lattice
                                                                              -- that would give 0 contribution anyway
                                                                              [Position i j | i <- [(-halfWidth) .. (halfWidth)], j <- [-(halfHeight) .. (halfHeight-1)]]
                                                                             )
                                                                             ++
                                                                             (map (\pos@(Position x y) -> ( positionToIndex pos
                                                                                , positionToIndex (Position x (y P.- 1))
                                                                                , 0 :+ (  plankConstant/2)
                                                                                )
                                                                                )
                                                                              -- removing min y from position range to keep operator within lattice
                                                                              -- that would give 0 contribution anyway
                                                                              [Position i j | i <- [(-halfWidth) .. (halfWidth)], j <- [-(halfHeight-1) .. halfHeight]]
                                                                             )
                                                                           )

        momSquaredOp = (momXOp `pow1p` 2) + (momYOp `pow1p` 2)
        posSquaredOp = (posXOp `pow1p` 2) + (posYOp `pow1p` 2)

        idOp :: SparseMatOperator (Complex Double) (CComplex CDouble) Position
        idOp = SparseMatOperator $! ES.fromList indexRange indexRange $! map (\pos@(Position _ _) -> (positionToIndex pos, positionToIndex pos, 1)) positionRange

        oscillationOperator :: Double -> SparseMatOperator (Complex Double) (CComplex CDouble) Position
        oscillationOperator dt =
              idOp + ( (0 :+ (0 P.- (dt P./ plankConstant))) .*
                       (
                         ( (((particleMass P./ 2) P.* (oscillatorFrequency P.^ 2)) :+ 0) .* posSquaredOp ) +
                         ( ((1 P./ (2 P.* particleMass)) :+ 0)                           .* momSquaredOp )
                       )
                     )

        initialModel :: SparseMatKet (Complex Double) (CComplex CDouble) Position
        initialModel = normalizeVector $! SparseMatKet $! ES.fromList indexRange 1 $! map ( \pos -> ( positionToIndex pos
                                                                                , 0
                                                                                , initialAmplitude pos
                                                                                )
                                                                                ) positionRange

        -- bitmap drawing
        drawing :: SparseMatKet (Complex Double) (CComplex CDouble) Position -> Picture
        drawing (SparseMatKet mat) =
            makePicture (2*halfWidth) (2*halfHeight)
                        1 1
                        ( \(x,y) -> let index = positionToIndex $! Position (round $! x P.* (P.fromIntegral halfWidth))
                                                                            (round $! y P.* (P.fromIntegral halfHeight))
                                        amp = mat ES.! (index, 0)
                                        den = double2Float $! normSquared amp
                                    in {-trace (show $ (index, amp, den)) $ -} rgb' den 0 0
                        )
            where
                normSquared cn = (magnitude cn) P.^ 2

        {-
        -- drawing with objects
        drawing :: SparseMatKet (Complex Double) (CComplex CDouble) Position -> Picture
        drawing (SparseMatKet mat) = {- trace "Drawing" $ -}
            pictures $ catMaybes $ map (\elem -> case elem of
                                           (index, 0, amp) -> let pos = indexToPosition index
                                                              in Just $! translate (fromIntegral $! posX pos) (fromIntegral $! posY pos) $!
                                                                 color (withAlpha (double2Float $! normSquared amp) red) $!
                                                                 circleSolid 1
                                           _ -> trace "OutOfBounds: " Nothing
                                       ) $! ES.toList $! mat
          where
                normSquared cn = (magnitude cn) P.^ 2
        -}

------------------------------------
-- simulate oscillation using sparse maps
------------------------------------

sparseModelSimulator :: OscillatorInABoxConfig
                     -> ( TimeStep -> SparseOperator (Complex Double) Position
                        , SparseKet (Complex Double) Position
                        , SparseKet (Complex Double) Position -> Picture
                        , OscillatorInABoxConfig
                        )
sparseModelSimulator config = (oscillationOperator, initialModel, drawing, config)
    where
        -- shorthands for config params
        positionRange = positionRangeConfig config
        halfWidth = halfWidthConfig config
        halfHeight = halfHeightConfig config

        plankConstant = plankConstantConfig config
        particleMass = particleMassConfig config
        oscillatorFrequency = oscillatorFrequencyConfig config

        initialAmplitude = initialAmplitudeConfig config

        -- position and momentum operators

        posXOp :: SparseOperator (Complex Double) Position
        posXOp = SparseOperator $! M.fromList $! map (\pos@(Position x _) -> ((pos,pos), ((fromIntegral x) :+ 0))) positionRange

        posYOp :: SparseOperator (Complex Double) Position
        posYOp = SparseOperator $! M.fromList $! map (\pos@(Position _ y) -> ((pos,pos), ((fromIntegral y) :+ 0))) positionRange

        momXOp :: SparseOperator (Complex Double) Position
        -- taking dx as 2
        -- <x'| P_x |a>  =  -ih ( d/dx(a)|_x' ) = - ( a(x'+1) - a(x'-1) ) * ih/2   = - (<x'+1|a> - <x'-1|a>) * ih / 2
        momXOp = SparseOperator $! M.fromList $! ( (map (\pos@(Position x y) -> ( (pos, Position (x P.+ 1) y), 0 :+ (- plankConstant/2) ))
                                                        -- removing max x from position range to keeo operator within lattice
                                                        -- that would give 0 contribution anyway
                                                        [Position i j | i <- [(-halfWidth) .. (halfWidth-1)], j <- [-(halfHeight) .. halfHeight]]
                                                   )
                                                   ++
                                                   (map (\pos@(Position x y) -> ( (pos, Position (x P.- 1) y), 0 :+ (  plankConstant/2) ))
                                                        -- removing min x from position range to keeo operator within lattice
                                                        -- that would give 0 contribution anyway
                                                        [Position i j | i <- [-(halfWidth-1) .. (halfWidth)], j <- [-(halfHeight) .. halfHeight]]
                                                   )
                                                 )

        momYOp :: SparseOperator (Complex Double) Position
        -- taking dy as 2
        -- <y'| P_y |a>  =  -ih ( d/dy(a)|_y' ) = - ( a(y'+1) - a(y'-1) ) * ih/2   = - (<y'+1|a> - <y'-1|a>) * ih / 2
        momYOp = SparseOperator $! M.fromList $! ( (map (\pos@(Position x y) -> ( (pos, Position x (y P.+ 1)), 0 :+ (- plankConstant/2) ))
                                                        -- removing max y from position range to keeo operator within lattice
                                                        -- that would give 0 contribution anyway
                                                        [Position i j | i <- [(-halfWidth) .. (halfWidth)], j <- [-(halfHeight) .. (halfHeight-1)]]
                                                   )
                                                   ++
                                                   (map (\pos@(Position x y) -> ( (pos, Position x (y P.- 1)), 0 :+ (  plankConstant/2) ))
                                                        -- removing min y from position range to keeo operator within lattice
                                                        -- that would give 0 contribution anyway
                                                        [Position i j | i <- [(-halfWidth) .. (halfWidth)], j <- [-(halfHeight-1) .. halfHeight]]
                                                   )

                                                 )

        momSquaredOp = (momXOp `pow1p` 2) + (momYOp `pow1p` 2)
        posSquaredOp = (posXOp `pow1p` 2) + (posYOp `pow1p` 2)

        idOp :: SparseOperator (Complex Double) Position
        idOp = SparseOperator $! M.fromList $! map (\pos@(Position _ _) -> ((pos,pos), 1)) positionRange

        oscillationOperator :: Double -> SparseOperator (Complex Double) Position
        oscillationOperator dt =
              idOp + ( (0 :+ (0 P.- (dt P./ plankConstant))) .*
                       (
                         ( (((particleMass P./ 2) P.* (oscillatorFrequency P.^ 2)) :+ 0) .* posSquaredOp ) +
                         ( ((1 P./ (2 P.* particleMass)) :+ 0)                           .* momSquaredOp )
                       )
                     )

        initialModel :: SparseKet (Complex Double) Position
        initialModel = normalizeVector $! SparseKet $! M.fromList $! map ( \pos -> (pos, initialAmplitude pos) ) positionRange

        -- bitmap drawing
        drawing :: SparseKet (Complex Double) Position -> Picture
        drawing (SparseKet coeffMap) =
            makePicture (2*halfWidth) (2*halfHeight)
                        1 1
                        ( \(x,y) -> let pos = Position (round $! x P.* (P.fromIntegral halfWidth))
                                                       (round $! y P.* (P.fromIntegral halfHeight))
                                        maybeAmp = M.lookup pos coeffMap
                                    in case maybeAmp of
                                           Nothing -> rgb' 0 0 0
                                           Just amp -> let den = double2Float $! normSquared amp
                                                       in {-trace (show $ (index, amp, den)) $ -} rgb' den 0 0
                        )
            where
                normSquared cn = (magnitude cn) P.^ 2

        {-
        -- drawing with objects
        drawing :: SparseKet (Complex Double) Position -> Picture
        drawing (SparseKet coeffMap) = {- trace "Drawing" $ -}
            pictures $ map (\((Position x y), amp) -> translate (fromIntegral x) (fromIntegral y) $!
                                                              color (withAlpha (double2Float $! normSquared amp) red) $!
                                                                  circleSolid 1
                           ) $! M.toList $! coeffMap
            where
                normSquared cn = (magnitude cn) P.^ 2
        -}

------------------------------------
-- Using repa arrays and matrices
-- Matrix multiplication turns out to be the bottleneck
-- Probably because it's not designed for sparse matrices
-- can work for around halfWidth=2, halfHeight=2
------------------------------------

arrayModelSimulator :: OscillatorInABoxConfig
                    -> ( TimeStep -> LatticeOperator (Complex Double) Position
                       , LatticeKet (Complex Double) Position
                       , LatticeKet (Complex Double) Position -> Picture
                       , OscillatorInABoxConfig
                       )
arrayModelSimulator config = (oscillationOperator, initialModel, drawing, config)
    where
        -- shorthands for config params
        positionRange = positionRangeConfig config
        halfWidth = halfWidthConfig config
        halfHeight = halfHeightConfig config

        indexRange = indexRangeConfig config
        positionToIndex = positionToIndexConfig config
        indexToPosition = indexToPositionConfig config

        plankConstant = plankConstantConfig config
        particleMass = particleMassConfig config
        oscillatorFrequency = oscillatorFrequencyConfig config

        initialAmplitude = initialAmplitudeConfig config

        -- position and momentum operators

        posXOp :: LatticeOperator (Complex Double) Position
        posXOp = LatticeOperator $! R.fromFunction (Z :. indexRange :. indexRange) (\(Z:.i:.j) -> if i==j
                                                                                then (fromIntegral $! posX $! indexToPosition i) :+ 0
                                                                                else 0
                                                                                )

        posYOp :: LatticeOperator (Complex Double) Position
        posYOp = LatticeOperator $! R.fromFunction (Z :. indexRange :. indexRange) (\(Z:.i:.j) -> if i==j
                                                                                then (fromIntegral $! posY $! indexToPosition i) :+ 0
                                                                                else 0
                                                                                )

        momXOp :: LatticeOperator (Complex Double) Position
        -- taking dx as 2
        -- <x'| P_x |a>  =  -ih ( d/dx(a)|_x' ) = - ( a(x'+1) - a(x'-1) ) * ih/2   = - (<x'+1|a> - <x'-1|a>) * ih / 2
        momXOp = LatticeOperator $! R.fromFunction (Z :. indexRange :. indexRange) (\(Z:.i:.j) -> let xo = posX $! indexToPosition j
                                                                                yo = posY $! indexToPosition j
                                                                                xn = posX $! indexToPosition i
                                                                                yn = posY $! indexToPosition i
                                                                                in
                                                                                case yn==yo of
                                                                                False -> 0
                                                                                True
                                                                                |xn == xo P.+ 1 -> 0 :+ (  plankConstant/2)
                                                                                |xn == xo P.- 1 -> 0 :+ (- plankConstant/2)
                                                                                |otherwise -> 0
                                                                                )

        momYOp :: LatticeOperator (Complex Double) Position
        -- taking dy as 2
        -- <y'| P_y |a>  =  -ih ( d/dy(a)|_y' ) = - ( a(y'+1) - a(y'-1) ) * ih/2   = - (<y'+1|a> - <y'-1|a>) * ih / 2
        momYOp = LatticeOperator $! R.fromFunction (Z :. indexRange :. indexRange) (\(Z:.i:.j) -> let xo = posX $! indexToPosition j
                                                                                yo = posY $! indexToPosition j
                                                                                xn = posX $! indexToPosition i
                                                                                yn = posY $! indexToPosition i
                                                                                in
                                                                                case xn==xo of
                                                                                False -> 0
                                                                                True
                                                                                |yn == yo P.+ 1 -> 0 :+ (  plankConstant/2)
                                                                                |yn == yo P.- 1 -> 0 :+ (- plankConstant/2)
                                                                                |otherwise -> 0
                                                                                )

        momSquaredOp = (momXOp `pow1p` 2) + (momYOp `pow1p` 2)
        posSquaredOp = (posXOp `pow1p` 2) + (posYOp `pow1p` 2)

        -- identity operator
        -- can make operator an instance of Unital, but don't know any one of writing the identity matrix in repa independent of dimensions
        idOp :: LatticeOperator (Complex Double) Position
        idOp = LatticeOperator $! R.fromFunction (Z :. indexRange :. indexRange) (\(Z:.i:.j) -> case i==j of
                                                                                False -> 0
                                                                                True -> 1
                                                                                )

        -- simulation operator
        oscillationOperator :: Double -> LatticeOperator (Complex Double) Position
        oscillationOperator dt =
              idOp + ( (0 :+ (0 P.- (dt P./ plankConstant))) .*
                       (
                         ( (((particleMass P./ 2) P.* (oscillatorFrequency P.^ 2)) :+ 0) .* posSquaredOp ) +
                         ( ((1 P./ (2 P.* particleMass)) :+ 0)                           .* momSquaredOp )
                       )
                     )

        -- initial ket
        initialModel :: LatticeKet (Complex Double) Position
        initialModel = LatticeKet $! R.fromFunction (Z :. indexRange) $! \(Z:.i) -> initialAmplitude $! indexToPosition i

        -- bitmap drawing
        drawing :: LatticeKet (Complex Double) Position -> Picture
        drawing (LatticeKet vector) =
            makePicture (2*halfWidth) (2*halfHeight)
                        1 1
                        ( \(x,y) -> let pos = Position (round $! x P.* (P.fromIntegral halfWidth))
                                                       (round $! y P.* (P.fromIntegral halfHeight))
                                        amp = vector R.! ( Z:.(positionToIndex pos) )
                                        den = double2Float $! normSquared amp
                                    in {-trace (show $ (pos, amp, den)) $ -} rgb' den 0 0
                        )
            where
                normSquared cn = (magnitude cn) P.^ 2

        {-
        -- drawing with objects
        drawing :: LatticeKet (Complex Double) Position -> Picture
        drawing (LatticeKet vector) = {- trace "Drawing" $ -}
          pictures $ map (\(Position x y) -> translate (fromIntegral x) (fromIntegral y) $!
                                    color (withAlpha (double2Float $ normSquared $! (vector R.! ( Z:.(positionToIndex (Position x y))))) red) $!
                                    circleSolid 1
                         ) $ positionRange
            where normSquared cn = (magnitude cn) P.^ 2
        -}

------------------------------------
-- simulate oscillation using index functions
-- doesn't work at the moment
-- no NormalizeVector instance for InfKet
------------------------------------

indexFunctionsSimulator :: OscillatorInABoxConfig
                        -> ( TimeStep -> InfOperator (Complex Double) Position
                           , InfKet (Complex Double) Position
                           , InfKet (Complex Double) Position -> Picture
                           , OscillatorInABoxConfig
                           )
indexFunctionsSimulator config = (oscillationOperator, initialModel, drawing, config)
    where
        -- shorthands for config params
        positionRange = positionRangeConfig config
        halfWidth = halfWidthConfig config
        halfHeight = halfHeightConfig config

        plankConstant = plankConstantConfig config
        particleMass = particleMassConfig config
        oscillatorFrequency = oscillatorFrequencyConfig config

        initialAmplitude = initialAmplitudeConfig config

        -- position and momentum operators
        posXOp :: InfOperator (Complex Double) Position
        posXOp = InfOperator $ \pos@(Position x _) -> M.fromList [(pos, (fromIntegral x) :+ 0) ]

        posYOp :: InfOperator (Complex Double) Position
        posYOp = InfOperator $ \pos@(Position _ y) -> M.fromList [(pos, (fromIntegral y) :+ 0)]

        momXOp :: InfOperator (Complex Double) Position
        -- taking dx as 2
        -- <x'| P_x |a>  =  -ih ( d/dx(a)|_x' ) = - ( a(x'+1) - a(x'-1) ) * ih/2   = - (<x'+1|a> - <x'-1|a>) * ih / 2
        momXOp = InfOperator $ \(Position x y) -> M.fromList [ ( Position (x P.- 1) y, (0 :+ (  plankConstant/2)) ),
                                                               ( Position (x P.+ 1) y, (0 :+ (- plankConstant/2)) )
                                                             ]

        momYOp :: InfOperator (Complex Double) Position
        -- taking dy as 2
        -- <y'| P_y |a>  =  -ih ( d/dy(a)|_y' ) = - ( a(y'+1) - a(y'-1) ) * ih/2   = - (<y'+1|a> - <y'-1|a>) * ih / 2
        momYOp = InfOperator $ \(Position x y) -> M.fromList [ ( Position x (y P.- 1), 0 :+ (  plankConstant/2) ),
                                                               ( Position x (y P.+ 1), 0 :+ (- plankConstant/2) )
                                                             ]

        momSquaredOp = (momXOp `pow1p` 2) + (momYOp `pow1p` 2)
        posSquaredOp = (posXOp `pow1p` 2) + (posYOp `pow1p` 2)

        scOp :: (Complex Double) -> InfOperator (Complex Double) Position
        scOp scalar = InfOperator $! \pos -> M.fromList [ (pos, scalar) ]

        oscillationOperator :: Double -> InfOperator (Complex Double) Position
        oscillationOperator dt = (scOp $! 1 :+ 0) + ((scOp $! 0 :+ (0 P.- (dt P./ plankConstant))) *
                                                     (
                                                       ( momSquaredOp * (scOp $ (1 P./ (2 P.* particleMass)) :+ 0) ) +
                                                       ( posSquaredOp * (scOp $ ((particleMass P./ 2) P.* oscillatorFrequency P.^ 2) :+ 0) )
                                                     )
                                                    )

        initialModel :: InfKet (Complex Double) Position
        initialModel = InfKet $! initialAmplitude

        -- bitmap drawing
        drawing :: InfKet (Complex Double) Position -> Picture
        drawing (InfKet coeffMap) = {- trace "Drawing" $ -}
            makePicture (2*halfWidth) (2*halfHeight)
                        1 1
                        ( \(x,y) -> let pos = Position (round $! x P.* (P.fromIntegral halfWidth))
                                                       (round $! y P.* (P.fromIntegral halfHeight))
                                        amp = coeffMap pos
                                        den = double2Float $! normSquared amp
                                    in {-trace (show $ (pos, amp, den)) $ -} rgb' den 0 0
                        )
            where
                normSquared cn = (magnitude cn) P.^ 2

        {-
        -- drawing with objects
        drawing :: InfKet (Complex Double) Position -> Picture
        drawing (InfKet coeffMap) = {- trace "Drawing" $ -}
            pictures $ map (\(Position x y) -> translate (fromIntegral x) (fromIntegral y) $!
                                               color (withAlpha (double2Float $! normSquared $! coeffMap (Position x y)) red) $!
                                               circleSolid 1
                           ) $ positionRange
            where normSquared cn = (magnitude cn) P.^ 2
        -}

------------------------------------------
-- Simulate
------------------------------------------

simulateOscillation :: (VectorSpace k v, SimpleOperator k v o, NormalizeVector v)
                    => (TimeStep -> o)
                    -> v
                    -> (v -> Picture)
                    -> OscillatorInABoxConfig
                    -> IO ()
simulateOscillation oscillationOperator initialModel drawing config = do
    -- return $! trace "Start Calculating Operator." ()
    -- keeping fps constant to avoid recalculation of operator
    operator <- return $! oscillationOperator $! timeFactor P./ (fromIntegral fps)
    -- return $! trace "Finish Calculating Operator." ()
    simulate (InWindow "Quantum Oscillation Simulator" (2*halfWidth, 2*halfHeight) (halfWidth, halfHeight))
             black
             fps
             initialModel
             drawing
             (simulation operator)

    where
        fps = fpsConfig config
        halfHeight = halfHeightConfig config
        halfWidth = halfWidthConfig config
        timeFactor = timeFactorConfig config

        simulation :: (VectorSpace k v, SimpleOperator k v o, NormalizeVector v)
                   => o -> ViewPort -> Float
                   -> v -> v
        simulation operator _ _ ket = {- trace "Simulating" $! -} normalizeVector $! operator `op` ket

textSimulateOscillation :: (VectorSpace k v, Show v, SimpleOperator k v o, NormalizeVector v)
                        => (TimeStep -> o)
                        -> v
                        -> OscillatorInABoxConfig
                        -> IO ()
textSimulateOscillation oscillationOperator initialModel config = do
    -- return $! trace "Start Calculating Operator." ()
    operator <- return $ oscillationOperator $! timeFactor P./ (fromIntegral fps)
    -- return $! trace "Finish Calculating Operator." ()
    go operator initialModel

    where
        fps = fpsConfig config
        timeFactor = timeFactorConfig config

        go :: ( VectorSpace k v, SimpleOperator k v o
              , Show v, NormalizeVector v
              )
            => o -> v -> IO ()
        go operator ket = do
                print ket
                threadDelay 1000000
                let nextMap = operator `op` ket
                go operator $! normalizeVector nextMap

--------------------------------
--Class Instances
--------------------------------

class NormalizeVector a where
    normalizeVector :: a -> a

instance NormalizeVector (SparseMatKet (Complex Double) (CComplex CDouble) Position) where
    normalizeVector ket@(SparseMatKet mat) = (1 / (norm :+ 0)) *^ ket
        where
            norm :: Double
            -- norm = realPart $! ES.norm mat
            -- normalizing with max of norms instead of sum of squares of norm to make the states more visible in simulation
            norm = ESD.maxCoeff $ ESD.convert normSquared $ ES.toMatrix mat

            normSquared cn = (magnitude cn) P.^ 2

instance NormalizeVector (SparseKet (Complex Double) Position) where
    normalizeVector ket@(SparseKet ketMap) = (1 / (norm :+ 0)) *^ ket
        where
            -- ketNormSquared = M.foldl' (P.+) 0 $! M.map (\coeff -> normSquared coeff) ketMap
            -- normalizing with max of norms instead of sum of squares of norm to make the states more visible in simulation
            maxNormSquared = M.foldl' (\a b -> max a b) 0 $! M.map (\coeff -> normSquared coeff) ketMap
            normSquared cn = (magnitude cn) P.^ 2
            norm :: Double
            norm = sqrt $! maxNormSquared

instance NormalizeVector (LatticeKet (Complex Double) Position) where
    normalizeVector ket@(LatticeKet array) = (1 / (norm :+ 0)) *^ ket
        where
            -- ketNormSquared = R.foldlAllS (P.+) 0 $! R.map (\coeff -> normSquared coeff) array
            -- normalizing with max of norms instead of sum of squares of norm to make the states more visible in simulation
            maxNormSquared = R.foldAllS (\a b -> max a b) 0 $! R.map (\coeff -> normSquared coeff) array
            normSquared cn = (magnitude cn) P.^ 2
            norm :: Double
            norm = sqrt $! maxNormSquared

{-
-- Will have to encode positionRange in Position to be able to traverse index map
instance NormalizeVector (InfKet (Complex Double) Position) where
    normalizeVector ket@(InfKet coeffMap) = (1 / (norm :+ 0)) *^ ket
        where
            -- ketNormSquared = foldl' (P.+) 0 $! map (normSquared .coeffMap) positionRange
            -- normalizing with max of norms instead of sum of squares of norm to make the states more visible in simulation
            maxNormSquared = foldl' (\a b -> max a b) 0 $! map (normSquared .coeffMap) positionRange
            normSquared cn = (magnitude cn) P.^ 2
            norm :: Double
            norm = sqrt $! maxNormSquared
-}

-- Hermitian Inner Product
instance (InnerProductSpace (Complex Double) (SparseKet (Complex Double) Position)) where
    (<.>) (SparseKet coeffA) (SparseKet coeffB) = M.foldl' (+) 0 $! M.intersectionWith (\ampA ampB -> ampA * (complexConjugate ampB)) coeffA coeffB

-------------------
-- For field instance of Complex Double

instance Euclidean (Complex Double) where degree _ = Just 1
instance Division (Complex Double) where (/) = (P./)
instance PID (Complex Double)
instance UFD (Complex Double)
instance GCDDomain (Complex Double)
instance IntegralDomain (Complex Double)
instance Unital (Complex Double) where one = 1 :+ 0
instance DecidableZero (Complex Double) where isZero = ((0 :+ 0) ==)
instance ZeroProductSemiring (Complex Double)
instance Commutative (Complex Double)
instance Multiplicative (Complex Double) where (*) = (P.*)
instance DecidableUnits (Complex Double) where recipUnit a = case a of
                                                (1 :+ 0) -> Just a
                                                (-1 :+ 0) -> Just a
                                                _ -> Nothing
instance UnitNormalForm (Complex Double)
instance Ring (Complex Double)
instance Monoidal (Complex Double) where zero = 0 :+ 0
instance Semiring (Complex Double)
instance Additive (Complex Double) where (+) = (P.+)
instance DecidableAssociates (Complex Double) where isAssociate = (==)
instance Rig (Complex Double)
instance Group (Complex Double) where (-) = (P.-)
instance LeftModule Natural (Complex Double) where
    n .* m = (fromIntegral $ toInteger n) P.* m
instance RightModule Natural (Complex Double) where
    m *. n = m P.* (fromIntegral n)
instance LeftModule Integer (Complex Double) where
    n .* m = (fromIntegral $ toInteger n) P.* m
instance RightModule Integer (Complex Double) where
    m *. n = m P.* (fromIntegral n)
instance Abelian (Complex Double)

--------------
-- For field instance of Double

instance Euclidean Double where degree _ = Just 1
instance Division Double where (/) = (P./)
instance PID Double
instance UFD Double
instance GCDDomain Double
instance IntegralDomain Double
instance Unital Double where one = 1
instance DecidableZero Double where isZero = (0 ==)
instance ZeroProductSemiring Double
instance Commutative Double
instance Multiplicative Double where (*) = (P.*)
instance DecidableUnits Double where recipUnit a = case a of
                                                       1 -> Just a
                                                       -1 -> Just a
                                                       _ -> Nothing
instance UnitNormalForm Double
instance Ring Double
instance Monoidal Double where zero = 0
instance Semiring Double
instance Additive Double where (+) = (P.+)
instance DecidableAssociates Double where isAssociate = (==)
instance Rig Double
instance Group Double where (-) = (P.-)
instance LeftModule Natural Double where
    n .* m = (fromIntegral $ toInteger n) P.* m
instance RightModule Natural Double where
    m *. n = m P.* (fromIntegral n)
instance LeftModule Integer Double where
    n .* m = (fromIntegral $ toInteger n) P.* m
instance RightModule Integer Double where
    m *. n = m P.* (fromIntegral n)
instance Abelian Double

---------------
-- For use by repa
instance Elt (Complex Double) where
    touch (a :+ b) = touch a >> touch b

max_stars_repo_name: sarthakbagaria/has-sci

max_issues_repo_name: sarthakbagaria/has-sci

max_forks_repo_name: sarthakbagaria/has-sci


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ScopedTypeVariables #-}

import Numeric.LinearAlgebra
import Control.Monad
import Control.Monad.Random
import Control.Applicative (liftA2)

data Weights = W { wBiases :: !(Vector Double) -- n
                 , wNodes :: !(Matrix Double)  -- n x m
                 }                             -- "m to n" layer

data Network :: * where
  O :: !Weights -> Network
  (:~) :: !Weights -> !Network -> Network

infixr 5 :~

-- | Generate random weights.

randomWeights :: MonadRandom m => Int -> Int -> m Weights
randomWeights i o = do
  seed1  :: Int <- getRandom
  seed2  :: Int <- getRandom
  let wB = randomVector seed1 Uniform o
      wN = uniformSample seed2 o (replicate i (-1, 1))
  return $ W wB wN

-- | Generate a random Network
randomNet :: MonadRandom m => Int -> [Int] -> Int -> m Network
randomNet i [] o = O <$> randomWeights i o
randomNet i (h:hs) o = liftA2 (:~) (randomWeights i h) (randomNet i hs o)

-- | Relu function
relu :: (Floating a, Ord a) => a -> a
relu = max 0

relu' :: (Floating a, Ord a) => a -> a
relu' x | x < 0 = 0
        | otherwise = 1

train :: Double -> Vector Double -> Vector Double -> Network -> Network
train = undefined

ext: hs

max_stars_repo_name: Lockne/neural-network-from-scratch

max_issues_repo_name: Lockne/neural-network-from-scratch

max_forks_repo_name: Lockne/neural-network-from-scratch


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE DefaultSignatures #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE PatternSynonyms #-}
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE UndecidableInstances #-}
-----------------------------------------------------------------------------
-- |
-- Copyright   :  (c) Edward Kmett 2010-2021
-- License     :  BSD3
-- Maintainer  :  ekmett@gmail.com
-- Stability   :  experimental
-- Portability :  GHC only
--
-----------------------------------------------------------------------------

module Numeric.AD.Mode
  (
  -- * AD modes
    Mode(..)
  , pattern KnownZero
  , pattern Auto
  ) where

import Numeric.Natural
import Data.Complex
import Data.Int
import Data.Ratio
import Data.Word

infixr 7 *^
infixl 7 ^*
infixr 7 ^/

class (Num t, Num (Scalar t)) => Mode t where
  type Scalar t
  type Scalar t = t

  -- | allowed to return False for items with a zero derivative, but we'll give more NaNs than strictly necessary
  isKnownConstant :: t -> Bool
  isKnownConstant _ = False

  asKnownConstant :: t -> Maybe (Scalar t)
  asKnownConstant _ = Nothing

  -- | allowed to return False for zero, but we give more NaN's than strictly necessary
  isKnownZero :: t -> Bool
  isKnownZero _ = False

  -- | Embed a constant
  auto  :: Scalar t -> t
  default auto :: (Scalar t ~ t) => Scalar t -> t
  auto = id

  -- | Scalar-vector multiplication
  (*^) :: Scalar t -> t -> t
  a *^ b = auto a * b

  -- | Vector-scalar multiplication
  (^*) :: t -> Scalar t -> t
  a ^* b = a * auto b

  -- | Scalar division
  (^/) :: Fractional (Scalar t) => t -> Scalar t -> t
  a ^/ b = a ^* recip b

  -- |
  -- @'zero' = 'lift' 0@
  zero :: t
  zero = auto 0

pattern KnownZero :: Mode s => s
pattern KnownZero <- (isKnownZero -> True) where
  KnownZero = zero

pattern Auto :: Mode s => Scalar s -> s
pattern Auto n <- (asKnownConstant -> Just n) where
  Auto n = auto n

instance Mode Double where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Float where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Int where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Integer where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Int8 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Int16 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Int32 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Int64 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Natural where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Word where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Word8 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Word16 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Word32 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Mode Word64 where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance RealFloat a => Mode (Complex a) where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

instance Integral a => Mode (Ratio a) where
  isKnownConstant _ = True
  asKnownConstant = Just
  isKnownZero x = 0 == x
  (^/) = (/)

max_stars_repo_name: msakai/ad

max_issues_repo_name: msakai/ad

max_forks_repo_name: msakai/ad


text:
import HVX
import Numeric.LinearAlgebra

-- declare a symbolic variable
x = EVar "x"
-- give x a value (a 4-element column vector)
vars = [("x", (4><1) [0,-3,1,2])]
-- define the expression to sudifferentiate: max(abs(x))
myexpr = hmax(habs(x))
-- compute the subgradient
mysubgrad = jacobianWrtVar myexpr vars "x"

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
module Duck.Statistics where

import Duck.Types

import Data.Ord (comparing)
import Data.List (maximumBy, sortBy)
import qualified Data.Vector.Unboxed as UVec
import Statistics.Regression
import Test.QuickCheck (generate)

-- |Generates a list of instance sizes to be benchmarked.
genSizeList :: (Int, Int) -> Int -> [Int]
genSizeList (lo, hi) iter = take iter [lo, lo + step .. hi]
  where step = (hi - lo) `div` iter

-- |Generates multiple instance of given sizes to be benchmarked.
genInstances :: (Cased a) => [Int] -> IO [a]
genInstances [] = return []
genInstances (x:xs) = (:) <$>
                      (generate $ genCase x) <*>
                      (genInstances xs)

-- |Packs and unpacks lists to perform a linear regression
-- on data and return a error value and a proportionality constant.
regress :: [Double] -> [Double] -> (Double, Double)
regress ana exp = (const UVec.! 0, conf)
  where (const, conf) = olsRegress [UVec.fromList (ana)] (UVec.fromList exp)

-- |A grouping function to be used with 'testGroup'. It
-- filters out the statistically unrelevant results.
relevant :: [(Double, a)] -> [(Double, a)]
relevant = filter ((>= 0.98) . fst)

-- |A grouping function to be used with 'testGroup'. It
-- keeps only the maximum confidence result.
maxConf :: [(Double, a)] -> [(Double, a)]
maxConf ls = [maximumBy (comparing fst) ls]

-- |A grouping function to be used with 'testGroup'. It
-- keeps only the top 'sz' maximum confidence results.
topConf :: Int -> [(Double, a)] -> [(Double, a)]
topConf sz ls = take sz $ reverse $ sortBy (comparing fst) ls

-- |A grouping function to be used with 'testGroup'. It
-- keeps all results.
allConf :: [(Double, a)] -> [(Double, a)]
allConf = id

max_stars_repo_name: gangsterveggies/duck-analyzer

max_issues_repo_name: gangsterveggies/duck-analyzer

max_forks_repo_name: gangsterveggies/duck-analyzer


text:
{-# LANGUAGE AllowAmbiguousTypes, FlexibleContexts, NamedFieldPuns #-}

module School.TestUtils
( CostFunc
, addClasses
, assertRight
, def
, diffCost
, diffInput
, doCost
, dummyList
, dummyMatrix
, empty
, fromLeft
, fromRight
, getRandDouble
, isSorted
, jTest
, matIndexes
, randomAffineParams
, randomMatrix
, randomMatrixL
, randomNNInts
, randomVector
, singleClassOutput
, testIOCatch
, testState
, unitCorrect
, weight1
, whenPrint
) where

import Conduit (ConduitM, liftIO, runConduit)
import Control.Exception (catch)
import Control.Monad (when)
import Data.Default.Class (def)
import Data.Either (either)
import Data.List (sort)
import Data.Void (Void)
import Numeric.LinearAlgebra ((><), (|>), Element, IndexOf, Matrix, R, Vector, accum, assoc, fromColumns, fromList, fromRows, size, sumElements, toColumns)
import School.Train.AppTrain (AppTrain, runAppTrain)
import School.Train.TrainState (TrainState)
import School.Types.Slinky (Slinky)
import School.Unit.CostFunction (CostFunction(..))
import School.Unit.CostParams (CostParams)
import School.Unit.Unit (Unit(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitParams (UnitParams(..))
import School.Unit.WeightDecay (weightDecay)
import School.Utils.Double (doubleRange)
import System.Exit (ExitCode(..))
import System.Random (getStdRandom, randomR)
import Test.QuickCheck.Monadic (PropertyM, assert)

testIOCatch :: IO a -> PropertyM IO ExitCode
testIOCatch action = liftIO $
  catch (action >> (return ExitSuccess))
        return

assertRight :: (b -> Bool)
            -> Either a b
            -> PropertyM IO ()
assertRight f = either (const $ assert False)
                       (\x -> assert $ f x)

testState :: ConduitM () Void (AppTrain R) b
          -> TrainState R
          -> PropertyM IO (Either String (b, TrainState R))
testState conduit state =
  liftIO . (runAppTrain state) . runConduit $ conduit

isSorted :: (Ord a) => [a] -> Bool
isSorted xs = xs == (reverse . sort $ xs)

dummyList :: (Num a) => Int -> Int -> [a]
dummyList r c = map fromIntegral [ r*i + j | i <- [1..r], j <- [1..c] ]

dummyMatrix :: Int -> Int -> Matrix Double
dummyMatrix r c = r >< c $ dummyList r c

getRandDouble :: IO Double
getRandDouble =
  getStdRandom (randomR doubleRange)

randomNNInts :: Int -> Int -> IO [Int]
randomNNInts maxInt n =
  sequence . (replicate n) $ getStdRandom (randomR (0, maxInt))

matMax :: Double
matMax = 1e3

randomMatrix :: Int -> Int -> IO (Matrix R)
randomMatrix nRows nCols = do
  let nEls = nRows * nCols
  inputList <- sequence . (replicate nEls) $ getStdRandom (randomR (-matMax, matMax))
  return $ (nRows >< nCols) inputList

randomMatrixL :: Double -> Int -> Int -> IO (Matrix R)
randomMatrixL limit nRows nCols = do
  let nEls = nRows * nCols
  inputList <- sequence . (replicate nEls) $ getStdRandom (randomR (-limit, limit))
  return $ (nRows >< nCols) inputList

randomVector :: Int -> IO (Vector R)
randomVector n = do
  inputList <- sequence . (replicate n) $ getStdRandom (randomR (-matMax, matMax))
  return $ n |> inputList

randomAffineParams :: Int
                   -> Int
                   -> IO (UnitParams R)
randomAffineParams fSize oSize = do
  affineBias <- randomVector oSize
  affineWeights <- randomMatrix oSize fSize
  return AffineParams { affineBias
                      , affineWeights
                      }

toPrint :: (Show a) => String -> a -> IO ()
toPrint tag val = putStrLn $ tag ++ " " ++ (show val)

whenPrint :: (Show a) => Bool -> [String] -> [a] -> IO ()
whenPrint cond tags vals = do
  when cond (sequence_ $ zipWith toPrint tags vals)

type CostFunc = UnitActivation R -> Double

jTest :: CostFunc
jTest (BatchActivation m) = sumElements m
jTest _ = 0

diffInput :: Unit R
          -> UnitParams R
          -> UnitActivation R
          -> Double
          -> IndexOf Matrix
          -> Double
diffInput unit params input eps idx =
  (jAdd - jSub) / (2*eps) where
    outAdd = apply unit params (alterInput eps idx input)
    outSub = apply unit params (alterInput (-eps) idx input)
    jAdd = jTest outAdd
    jSub = jTest outSub

fromRight :: b -> Either a b -> b
fromRight b = either (const b) id

fromLeft :: a -> Either a b -> a
fromLeft b = either id (const b)

diffCost :: CostFunction R m
         -> UnitActivation R
         -> Double
         -> IndexOf Matrix
         -> Slinky CostParams
         -> Double
diffCost costFunc input eps idx costParams = let
  jAdd = computeCost costFunc (alterInput eps idx input) costParams
  jSub = computeCost costFunc (alterInput (-eps) idx input) costParams
  in (fromRight 0 jAdd - fromRight 0 jSub) / (2*eps)

alterInput :: AlterInput
alterInput change idx (BatchActivation m) =
  BatchActivation $ accum m (+) [(idx, change)]
alterInput _ _ _ = ApplyFail "alterInput error"

type AlterInput = Double
               -> IndexOf Matrix
               -> UnitActivation R
               -> UnitActivation R

matIndexes :: Int -> Int -> [IndexOf Matrix]
matIndexes r c = [ (j, k) | j <- [0..r-1], k <- [0..c-1] ]

empty :: (Element a) => Matrix a
empty = (0><0) []

weight1 :: (Monad m) => CostFunction R m
weight1 = weightDecay 1

doCost :: (Element a, Num a)
       => CostFunction a m
       -> UnitActivation a
       -> Slinky CostParams
       -> (a, UnitGradient a)
doCost costFunction activation params =
  fromRight (0, BatchGradient empty) result where
    result = do
      cost <- computeCost costFunction activation params
      grad <- derivCost costFunction activation params
      return (cost, grad)

appCol :: (Element a)
       => Vector a
       -> Matrix a
       -> Matrix a
appCol v =
  fromColumns . (flip (++) $ [v]) . toColumns

addClasses :: (Element a, Num a)
           => [Int]
           -> Matrix a
           -> Matrix a
addClasses classes mat = appCol v mat where
  v = fromList . (map fromIntegral) $ classes

singleClassOutput :: Int -> [Int] -> UnitActivation R
singleClassOutput nClasses =
  BatchActivation . fromRows . (map (\idx -> assoc nClasses 0 [(idx, 1)]))

unitCorrect :: [Int] -> Unit R
unitCorrect classes = Unit { apply, deriv } where
  apply _ (BatchActivation input) =
    singleClassOutput n classes
    where (_, n) = size input
  apply _ _ = ApplyFail ""
  deriv = undefined

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
module STCR2Z2T0S0Edge where

import           Control.Arrow
import           Control.Monad
import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.Edge
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Types

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:thetaFreqsStr:scaleFreqsStr:histFilePath:numIterationStr:writeSourceFlagStr:edgeFilePath:numNoisePointStr:scaleFactorStr:useFFTWWisdomFlagStr:fftwWisdomFileName: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
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Double
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numNoisePoint = read numNoisePointStr :: Int
      scaleFactor = read scaleFactorStr :: Double
      useFFTWWisdomFlag = read useFFTWWisdomFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z2T0S0Edge"
      fftwWisdomFilePath = folderPath </> fftwWisdomFileName
  createDirectoryIfMissing True folderPath
  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
          thetaFreqs
          thetaFreqs
          scaleFreqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scaleFreqs
             , L.length thetaFreqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z2T0S0 <-
    computeUnboxedP $
    computeR2Z2T0S0ArrayRadial
      (PinwheelHollow0 4)
      radialArr
      numPoint
      numPoint
      1
      maxScale
      thetaFreqs
      scaleFreqs
      thetaFreqs
      scaleFreqs
  plan <-
    makeR2Z2T0S0Plan emptyPlan useFFTWWisdomFlag fftwWisdomFilePath arrR2Z2T0S0
  xs <- parseEdgeFile edgeFilePath
  randomPonintSet <- generateRandomPointSet numNoisePoint numPoint numPoint
  let (centerX, centerY) =
        join (***) (\x -> div x . L.length $ xs) .
        L.foldl' (\(a, b) (R2S1RPPoint (c, d, _, _)) -> (a + c, b + d)) (0, 0) $
        xs
      ys =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) ->
             (R2S1RPPoint
                ( round $ (fromIntegral $ a - centerX) / scaleFactor
                , round $ (fromIntegral $ b - centerX) / scaleFactor
                , c
                , d)))
          xs
      zs =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) ->
             (R2S1RPPoint (a - center numPoint, b - center numPoint, c, d)))
          randomPonintSet
      points = ys L.++ zs
  let bias = computeBiasR2T0S0 numPoint numPoint thetaFreqs scaleFreqs points
      eigenVec =
        computeInitialEigenVectorR2T0S0
          numPoint
          numPoint
          thetaFreqs
          scaleFreqs
          thetaFreqs
          scaleFreqs
          points
  powerMethodR2Z2T0S0
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    numScale
    scaleFreqs
    maxScale
    arrR2Z2T0S0
    numIteration
    writeSourceFlag
    (printf
       "_%d_%d_%.2f_%.2f"
       (round maxScale :: Int)
       (round tao :: Int)
       thetaSigma
       scaleSigma)
    bias
    eigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
-----------------------------------------------------------------------------
-- |
-- Module     : Numeric.LinearAlgebra.Matrix.Statistics
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Basic multivariate statistics.
--

module Numeric.LinearAlgebra.Matrix.Statistics (
    defaultCovUplo,

    -- * Immutable interface
    cov,
    covWithMean,
    weightedCov,
    weightedCovWithMean,

    -- * Mutable interface
    covTo,
    covWithMeanTo,
    weightedCovTo,
    weightedCovWithMeanTo,

    ) where

import Control.Monad( when )
import Control.Monad.ST( ST )
import Data.List( foldl' )
import Text.Printf( printf )

import Numeric.LinearAlgebra.Types

import Numeric.LinearAlgebra.Vector( Vector, RVector )
import qualified Numeric.LinearAlgebra.Vector as V

import Numeric.LinearAlgebra.Matrix.Base( Matrix )
import Numeric.LinearAlgebra.Matrix.STBase( STMatrix )
import qualified Numeric.LinearAlgebra.Matrix.STBase as M
import qualified Numeric.LinearAlgebra.Matrix.Herm as M

-- | Returns the default storage scheme for covariance matrices.
defaultCovUplo :: Uplo
defaultCovUplo = Lower

-- | Returns the sample covariance matrix as a hermitian matrix with storage
-- scheme equal to 'defaultCovUplo'.  The first argument gives the dimension
-- of the vectors.
cov :: (BLAS3 e)
          => Int -> CovMethod -> [Vector e] -> Herm Matrix e
cov p t xs = M.hermCreate $ do
    c <- Herm uplo `fmap` M.new_ (p,p)
    covTo c t xs
    return c
  where
    uplo = defaultCovUplo

-- | Given the pre-computed mean, returns the sample covariance matrix
-- with storage scheme equal to 'defaultCovUplo'.
covWithMean :: (BLAS3 e)
                  => Vector e -> CovMethod -> [Vector e] -> Herm Matrix e
covWithMean mu t xs = M.hermCreate $ do
    p <- V.getDim mu
    c <- Herm uplo `fmap` M.new_ (p,p)
    covWithMeanTo c mu t xs
    return c
  where
    uplo = defaultCovUplo

-- | Returns the weighed sample covariance matrix with storage scheme equal
-- to 'defaultCovUplo'. The first argument gives the dimension of the vectors.
weightedCov :: (BLAS3 e)
                  => Int -> CovMethod -> [(Double, Vector e)] -> Herm Matrix e
weightedCov p t wxs = M.hermCreate $ do
    c <- Herm uplo `fmap` M.new_ (p,p)
    weightedCovTo c t wxs
    return c
  where
    uplo = defaultCovUplo

-- | Given the pre-computed mean, returns the weighed sample covariance matrix
-- with storage scheme equal to 'defaultCovUplo'.
weightedCovWithMean :: (BLAS3 e)
                    => Vector e -> CovMethod -> [(Double, Vector e)]
                    -> Herm Matrix e
weightedCovWithMean mu t wxs = M.hermCreate $ do
    p <- V.getDim mu
    c <- Herm uplo `fmap` M.new_ (p,p)
    weightedCovWithMeanTo c mu t wxs
    return c
  where
    uplo = defaultCovUplo

-- | Computes and copies the sample covariance matrix to the given
-- destination.
covTo :: (RVector v, BLAS3 e)
      => Herm (STMatrix s) e -> CovMethod -> [v e] -> ST s ()
covTo c@(Herm _ a) t xs = do
    (p,_) <- M.getDim a
    mu <- V.new p 1
    V.meanTo mu xs
    covWithMeanTo c mu t xs

-- | Given the pre-computed mean, computes and copies the sample covariance
-- matrix to the given destination.
covWithMeanTo :: (RVector v1, RVector v2, BLAS3 e)
              => Herm (STMatrix s) e -> v1 e -> CovMethod -> [v2 e] -> ST s ()
covWithMeanTo c@(Herm _ a) mu t xs = do
    (ma,na) <- M.getDim a
    p <- V.getDim mu

    when ((ma,na) /= (p,p)) $ error $
        printf ("covWithMeanTo"
                ++ " (Herm _ <matrix with dim (%d,%d)>)"
                ++ " <vector with dim %d>"
                ++ " _ _"
                ++ ": dimension mismatch")
               ma na p

    xt <- M.new_ (p,n)
    M.withColsM xt $ \xs' ->
        sequence_ [ V.subTo x' mu x
                  | (x,x') <- zip xs xs'
                  ]
    M.hermRankKUpdateM_ (1/df) NoTrans xt 0 c
  where

    n = length xs
    df = fromIntegral $ case t of { MLCov -> n ; UnbiasedCov -> n - 1 }

-- | Computes and copies the weighed sample covariance matrix to the
-- given destination.
weightedCovTo :: (RVector v, BLAS3 e)
              => Herm (STMatrix s) e -> CovMethod -> [(Double, v e)] -> ST s ()
weightedCovTo c@(Herm _ a) t wxs = do
    (p,_) <- M.getDim a
    mu <- V.new p 1
    V.weightedMeanTo mu wxs
    weightedCovWithMeanTo c mu t wxs

-- | Given the pre-computed mean, computes and copies the weighed sample
-- covariance matrix to the given destination.
weightedCovWithMeanTo :: (RVector v1, RVector v2, BLAS3 e)
                      => Herm (STMatrix s) e
                      -> v1 e -> CovMethod -> [(Double, v2 e)]
                      -> ST s ()
weightedCovWithMeanTo c@(Herm _ a) mu t wxs = do
    (ma,na) <- M.getDim a
    p <- V.getDim mu

    when ((ma,na) /= (p,p)) $ error $
        printf ("weightedCovWithMeanTo"
                ++ " (Herm _ <matrix with dim (%d,%d)>):"
                ++ " <vector with dim %d>"
                ++ " _ _"
                ++ " dimension mismatch")
               ma na p

    xt <- M.new_ (p,n)
    M.withColsM xt $ \xs' ->
        sequence_ [  V.subTo x' mu x
                  >> V.scaleM_ (realToFrac $ sqrt (w / invscale)) x'
                  |  (w,x,x') <- zip3 ws xs xs'
                  ]
    M.hermRankKUpdateM_ 1 NoTrans xt 0 c
  where
    (ws0,xs) = unzip wxs
    w_sum = foldl' (+) 0 ws0
    ws = if w_sum == 0 then ws0 else map (/w_sum) ws0
    w2s_sum = foldl' (+) 0 $ map (^^(2::Int)) ws
    invscale = case t of
                   MLCov -> 1
                   UnbiasedCov -> (1 - w2s_sum)
    n = length ws0

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
{-# LANGUAGE GADTs #-}
module Math.IRT.Internal.LogLikelihood where

import Numeric.AD (Mode, Scalar)
import Statistics.Distribution (Distribution)

class (Distribution d) => LogLikelihood d where
    logLikelihood :: (Mode a, Floating a, Scalar a ~ Double) => Bool -> d -> a -> a

logLikeFunc :: (Distribution d, Mode a, Floating a, Scalar a ~ Double) => (d -> a -> a) -> Bool -> d -> a -> a
logLikeFunc f True  = (log .) . f
logLikeFunc f False = ((log . (1-)) .) . f

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module LispVals (
    Frame,
    Env,
    ThrowsError,
    EnvThrowsError,
    IOThrowsError,
    LispVal(..),
    LispError(..),
    LispEval
  ) where

import Data.Ratio (Ratio, numerator, denominator)
import Data.Complex (Complex((:+)))
import qualified Data.Vector as V (toList, Vector)
import qualified Data.Map.Strict as M (Map)
import Text.ParserCombinators.Parsec (ParseError)
import Control.Monad.Error (ErrorT, Error, noMsg, strMsg, MonadError, throwError, catchError)
import Control.Monad.State (State, StateT)
import Control.Monad.Cont (ContT(..), runContT)
import Control.Monad.Trans (lift)
import System.IO (Handle)
import Data.IORef (IORef)

type Frame = M.Map String LispVal
type Env = [IORef Frame]
type ThrowsError = Either LispError
-- ErrorT e m a ~ m (Either e a)
-- State s a ~ s -> (a, s)
-- ErrorT e (State s a) ~ s -> (Either e a, s)
type EnvThrowsError = ErrorT LispError (State Env)
-- StateT s m a ~ s -> m (a, s)
-- ErrorT e (StateT s m) a ~ s -> m (Either e a, s)
-- StateT s (ErrorT e m) a ~ s -> m (Either e (a, s))
type IOThrowsError = ErrorT LispError (StateT Env IO)

type LispEval = ContT LispVal IOThrowsError LispVal
-- mtl doesn't define a MonadError instances for ContT
-- http://stackoverflow.com/questions/10742151/why-cant-contt-be-made-an-instance-of-monaderror
instance MonadError e m => MonadError e (ContT r m) where
  throwError = lift . throwError
  catchError op h = ContT $ \k -> catchError (runContT op k) (\e -> runContT (h e) k)

data LispVal = Void
        | Atom String
        | List [LispVal]
        | Vector (V.Vector LispVal)
        | DottedList [LispVal] LispVal
        | Number Integer
        | String String
        | Bool Bool
        | Char Char
        | Float Float
        | Complex (Complex Float)
        | Ratio (Ratio Integer)
        | Port Handle
        | Continuation (LispVal -> LispEval)
        | PrimitiveFunc ([LispVal] -> ThrowsError LispVal)
        | IOFunc ([LispVal] -> IOThrowsError LispVal)
        | SpecialFormFunc ([LispVal] -> LispEval)
        | Func
            { params :: [String]
            , vararg :: Maybe String
            , body :: [LispVal], closure :: Env}

instance Show LispVal where show = showVal
showVal :: LispVal -> String
showVal Void = "<void>"
showVal (Continuation _) = "<continuation>"
showVal (String contents) = "\"" ++ unescapeString contents ++ "\""
showVal (Atom name) = name
showVal (Number contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (Char contents) = "#\\" ++ case contents of
                            '\n' -> "newline"
                            ' ' -> "space"
                            _ -> [contents]
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (DottedList h t) = "(" ++ unwordsList h ++ " . " ++ showVal t ++ ")"
showVal (Vector contents) =
    "#(" ++ (unwordsList . V.toList) contents ++ ")"
showVal (Float contents) = show contents
showVal (Complex (r :+ i))
    | i > 0 = show r ++ "+" ++ show i ++ "i"
    | i < 0 = show r ++ show i ++ "i"
    | otherwise = show r
showVal (Ratio contents) = show (numerator contents) ++ "/" ++ show (denominator contents)
showVal (Port _) = "<port>"
showVal (PrimitiveFunc _) = "<primitive>"
showVal (SpecialFormFunc _) = "<special form primitive>"
showVal (IOFunc _) = "<io primitive>"
showVal (Func {params = args, vararg = varargs}) =
  "(lambda (" ++ unwords (map show args) ++
     (case varargs of
        Nothing -> ""
        Just arg -> " . " ++ arg) ++ ") ...)"

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

showError :: LispError -> String
showError (UnboundVar message varname) = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func) = message ++ ": " ++ show func
showError (NumArgs expected found) = "Expected " ++ show expected
                                  ++ " args; found values " ++ unwordsList found
showError (TypeMismatch expected found) = "Invalid type: expected " ++ expected
                                       ++ ", found " ++ show found
showError (Parser parseErr) = "Parse error at " ++ show parseErr
showError (Unspecified message) = "Unspecified operation: " ++ message
showError (WrongClause) = "No matching clause"
showError (Default message) = "Unknown error: " ++ message

instance Show LispError where show = showError

instance Error LispError where
     noMsg = Default "An error has occurred"
     strMsg = Default

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

unescapeString :: String -> String
unescapeString [] = []
unescapeString (x:xs)
    | x == '\n' = "\\n" ++ unescapeString xs
    | x == '\t' = "\\t" ++ unescapeString xs
    | x == '\r' = "\\r" ++ unescapeString xs
    | x == '\"' = "\\\"" ++ unescapeString xs
    | x == '\\' = "\\\\" ++ unescapeString xs
    | otherwise = x : unescapeString xs

max_stars_repo_name: tismith/tlisp

max_issues_repo_name: tismith/tlisp

max_forks_repo_name: tismith/tlisp


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}

{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise #-}

-- The unused variables are left as documentation
{-# OPTIONS -Wno-unused-local-binds #-}

{-|
Module      : ControlSystems.DynamicSystems.Conversion
Description : Converts continuous time systems to discrete time, using static vectors.
Copyright   : (c) Ryan Orendorff, 2020
License     : BSD3
Stability   : experimental
-}
module ControlSystems.DynamicSystems.Conversions.Static (
  c2d
) where

import GHC.TypeLits
import Numeric.LinearAlgebra.Static

-- | Convert a continuous time linear differential system to a discrete form,
-- using state space models, using type level naturals.
--
-- Specifically note this implementation against
-- `ControlSystems.DynamicSystems.Conversions.c2d`, which takes many more
-- code to handle the complexities of creating a block matrix of the right
-- size and then splitting it apart again.
c2d :: forall m n. (KnownNat m, KnownNat n) =>
       Double         -- ^ The sampling time
    -> L n n          -- ^ The state transition matrix, often called A
    -> L n m          -- ^ The state input matrix, often called B
    -> (L n n, L n m) -- ^ The discretized A & B system
c2d t a b = discreteSystem
  where
    block = a ||| b  -- [ A B ]
              ===    -- [ 0 0 ]
               0

    -- Calculate the discrete time block matrix, which is the following
    -- ⌈ A_d B_d ⌉
    -- ⌊  0   I  ⌋
    exp_block = expm (konst t * block)

    -- Finally we can extract the correct submatrices assuming we pull out
    -- the right pieces. Note that we do not need to pass the indices here
    -- directly; the type system ensures that we pull out the correct A_d
    -- and B_d parts of the exp_block matrix.
    discreteSystem@(a_d, b_d) = splitCols $ fst (splitRows exp_block)

max_stars_repo_name: ryanorendorff/convex

max_issues_repo_name: ryanorendorff/convex

max_forks_repo_name: ryanorendorff/convex


text:
{-# LANGUAGE GADTs, Rank2Types, DataKinds, TypeOperators #-}

module Lib.Utility.Mnist.IO
(
    MnistTrainingSize, MnistTestSize, MnistData, TrainingData,
    getMnistData
) where

import Network.HTTP as HTTP
import Data.ByteString.Lazy as Lazy (ByteString, readFile, writeFile, empty)
import Data.ByteString.Lazy.Char8 as Char8 (pack)
import Codec.Compression.GZip as GZip
import qualified Data.IDX as IDX
import Control.Monad.Trans.Maybe
import Data.Vector.Unboxed (Vector, toList)
import qualified Numeric.LinearAlgebra as HMatrix
import qualified Numeric.LinearAlgebra.Static as Static
import System.IO.Error
import qualified Data.Serialize as Serialize
import Control.Arrow
import Control.Applicative
import GHC.TypeLits
import GHC.Generics (Generic)
import Data.Proxy

type TrainingData n = (Static.L n 784, Static.L n 10)

type MnistTrainingSize = 60000
type MnistTestSize = 10000
type MnistData = (TrainingData MnistTrainingSize, TrainingData MnistTestSize)

type SerializedMnistData = ([(Int, [Int])], [(Int, [Int])])
type LabeledIntData = [(Int, Vector Int)]

getMnistData :: MaybeT IO MnistData
getMnistData = do
    loadResult <- loadMnistData
    (tr, te) <- case loadResult of
        Just x -> return x
        Nothing -> do
            () <- MaybeT (Just <$> putStrLn "Try fetch for load fail")
            x <- fetchMnistData
            () <- saveMnistData x
            MaybeT $ return $ parseMatrix $ fromLabeledIntData2Matrix x
    return (tr, te)

loadMnistData :: MaybeT IO (Maybe MnistData)
loadMnistData = MaybeT $ do
    putStrLn $ "Loading from " ++ fpath
    matrixMaybe <- HMatrix.loadMatrix' fpath -- 70000><785
    let loadResult = matrixMaybe >>= parseMatrix
    return $ Just loadResult

-- toOneHot :: HMatrix.Vector Double -> HMatrix.Matrix Double
toOneHot :: forall n. KnownNat n => Static.L n 1 -> Static.L n 10
toOneHot label = Static.build f
    where m0 = Static.extract $ (Static.konst 0 :: Static.L 10 10)
          f r c = if c == (atIndex r 0) then 1.0 else 0
          atIndex r c = Static.unwrap label `HMatrix.atIndex` (truncate r, truncate c)

parseMatrix :: HMatrix.Matrix Double -> Maybe MnistData
parseMatrix matrix = do
        matrix' <- Static.create matrix :: Maybe (Static.L (MnistTrainingSize + MnistTestSize) 785)
        let (label, image) = Static.splitCols matrix'
        let (trL, teL) = Static.splitRows label :: (Static.L MnistTrainingSize 1, Static.L MnistTestSize 1)
        let (trD, teD) = Static.splitRows $ image / 255
        return ((trD, toOneHot trL), (teD, toOneHot teL))

fromLabeledIntData2Matrix :: (LabeledIntData, LabeledIntData) -> HMatrix.Matrix Double
fromLabeledIntData2Matrix (tr, te) = matrix
    where transform = map (second toList)
          tr' = transform tr
          te' = transform te
          join = map (uncurry (:))
          matrix = HMatrix.fromLists $ map (map realToFrac) $ join tr' ++ join te'

saveMnistData :: (LabeledIntData, LabeledIntData) -> MaybeT IO ()
saveMnistData labeledIntData = MaybeT $ do
    putStrLn $ "Saving to " ++ fpath
    let matrix = fromLabeledIntData2Matrix labeledIntData
    let saveAction = HMatrix.saveMatrix fpath "%lf" matrix
    () <- catchIOError saveAction (\_ -> do putStrLn "SaveError"; return ())
    return $ Just ()

fetchMnistData :: MaybeT IO (LabeledIntData, LabeledIntData)
fetchMnistData = do
    (Data_ trD) <- fetchIDX trainImagesPath decodeIDXData
    (Labels_ trL) <- fetchIDX trainLabelsPath decodeIDXLabels
    (Data_ teD) <- fetchIDX testImagesPath decodeIDXData
    (Labels_ teL) <- fetchIDX testLabelsPath decodeIDXLabels
    MaybeT $ return $ do
        tr <- IDX.labeledIntData trL trD
        te <- IDX.labeledIntData teL teD
        return (tr, te)

data IDXW = Data_ IDX.IDXData | Labels_ IDX.IDXLabels

decodeIDXData :: ByteString -> Maybe IDXW
decodeIDXData = fmap Data_ . IDX.decodeIDX
decodeIDXLabels :: ByteString -> Maybe IDXW
decodeIDXLabels = fmap Labels_ . IDX.decodeIDXLabels
fetchIDX :: String -> (ByteString -> Maybe IDXW) -> MaybeT IO IDXW
fetchIDX path decodeF = MaybeT $ do
    let url = toUrl path
    putStrLn $ "Downloading from " ++ url
    idxBytes <- getFileFromUrl url
    let idx = decodeF idxBytes
    return idx

baseUrl = "http://yann.lecun.com/exdb/mnist/"
trainImagesPath = "train-images-idx3-ubyte.gz"
trainLabelsPath = "train-labels-idx1-ubyte.gz"
testImagesPath = "t10k-images-idx3-ubyte.gz"
testLabelsPath = "t10k-labels-idx1-ubyte.gz"
fpath = "mnist.pickle"
toUrl path = baseUrl ++ path

getFileFromUrl :: String -> IO ByteString
getFileFromUrl url = do
    response <- simpleHTTP (getRequest url)
    body <- getResponseBody response
    let bytes = Char8.pack body
    return $ GZip.decompress bytes


text:
{-# LANGUAGE CPP  #-}
#if defined(__GLASGOW_HASKELL__) && __GLASGOW_HASKELL__ >= 702
{-# LANGUAGE Safe #-}
#endif
-------------------------------------------------------------------------------
-- |
-- Module      : OGLS.Engine.Math.Vectors
-- description : 3D and 4D vectors arithmetic
-- Copyright   : (c) Adam, 2017
-- License     : MIT
-- Maintainer  : awkure@protonmail.ch
-- Stability   : unstable
-- Portability : POSIX
-------------------------------------------------------------------------------

module OGLS.Engine.Math.Instances where

import Data.Hashable              ( Hashable (..)         )
import Control.Monad.Fix          ( MonadFix (..)         )
import Control.Monad.Zip          ( MonadZip (..)         )
import Control.Applicative        ( liftA2                )
import Data.Functor.Bind          ( Apply (..), Bind (..) )
import Data.Semigroup.Foldable    ( Foldable1 (..)        )
import Data.Semigroup.Traversable ( Traversable1 (..)     )
import Data.Complex               ( Complex (..)          )
import Data.Semigroup             ( (<>)                  )
import Data.HashMap.Lazy as HashMap

instance (Hashable k, Eq k) => Apply (HashMap k) where
  (<.>) = HashMap.intersectionWith id

instance (Hashable k, Eq k) => Bind (HashMap k) where
  m >>- f = HashMap.fromList $ do
    (k, a) <- HashMap.toList m
    case HashMap.lookup k (f a) of
      Just b -> [(k,b)]
      Nothing -> []

instance Apply Complex where
  (a :+ b) <.> (c :+ d) = a c :+ b d

instance Bind Complex where
  (a :+ b) >>- f = a' :+ b' where
    a' :+ _  = f a
    _  :+ b' = f b
  {-# INLINE (>>-) #-}

instance MonadZip Complex where
  mzipWith = liftA2

instance MonadFix Complex where
  mfix f = (let a :+ _ = f a in a) :+ (let _ :+ a = f a in a)

instance Foldable1 Complex where
  foldMap1 f (a :+ b) = f a <> f b
  {-# INLINE foldMap1 #-}

instance Traversable1 Complex where
  traverse1 f (a :+ b) = (:+) <$> f a <.> f b
  {-# INLINE traverse1 #-}

max_stars_repo_name: awkure/ogls

max_issues_repo_name: awkure/ogls

max_forks_repo_name: awkure/ogls


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE UndecidableInstances #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Graphics.Rendering.Plot.Figure
-- Copyright   :  (c) A. V. H. McPhail 2010
-- License     :  BSD3
--
-- Maintainer  :  haskell.vivian.mcphail <at> gmail <dot> com
-- Stability   :  provisional
-- Portability :  portable
--
-- Creation and manipulation of 'Figure's
--
-- The same problem of leaked instances as at <http://hackage.haskell.org/packages/archive/graphviz/2999.10.0.1/doc/html/Data-GraphViz-Commands.html#t%3AGraphvizCanvas> occurs here.
--
--
-- /with/, /set/, /clear/, /new/, and /add/ are the operations that can
-- be performed on various elements of a figure.
--
-- /glib/\//data-accessor/ abstractions (verbs/modifiers) are planned for future implementations

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

module Graphics.Rendering.Plot.Figure (
                                       module Data.Colour.Names
                                       -- * Top level operation
                                      , Figure(), FigureState()
                                       -- * Default options
                                      , withTextDefaults
                                      , withLineDefaults
                                      , withPointDefaults
                                      , withBarDefaults
                                      -- * Figures
                                      , newFigure
                                      -- ** Formatting
                                      , setBackgroundColour
                                      , setFigurePadding
                                      , withTitle
                                      , withSubTitle
                                      , setPlots
                                      , withPlot, withPlots
                                      -- * Sub-plots
                                      ,  Plot()
                                      -- ** Colour
                                      , setPlotBackgroundColour
                                      -- ** Plot elements
                                      , Border
                                      , setBorder
                                      , setPlotPadding
                                      , withHeading
                                      -- ** Series data
                                      , Function(), VectorFunction(), Series(), MinMaxSeries(), ErrorSeries()
                                      , Surface()
                                      , SeriesLabel()
                                      , Abscissa(), Ordinate(), Dataset()
                                      , FormattedSeries(), SeriesType(..)
                                      , line, point, linepoint
                                      , impulse, step
                                      , area
                                      , bar
                                      , hist
                                      , candle, whisker
                                      , setDataset
                                      -- * Annotations
                                      , Location, Head, Fill
                                      , Annote()
                                      , arrow
                                      , oval
                                      , rect
                                      , glyph
                                      , text
                                      , cairo
                                      , withAnnotations
                                      -- ** Plot type
                                      , setSeriesType
                                      , setAllSeriesTypes
                                      -- ** Formatting
                                      , PlotFormats()
                                      , withSeriesFormat
                                      , withAllSeriesFormats
                                      -- * Range
                                      , Scale(..)
                                      , setRange
                                      , setRangeFromData
                                      -- * Axes
                                      , Axis
                                      , AxisType(..),AxisSide(..),AxisPosn(..)
                                      , clearAxes
                                      , clearAxis
                                      , addAxis
                                      , withAxis
                                      -- * BarSetting
                                      , BarSetting(..)
                                      , barSetting
                                      -- * Data Sampling
                                      , SampleData
                                      , sampleData
                                      -- * Legend
                                      , Legend
                                      , LegendBorder
                                      , LegendLocation(..), LegendOrientation(..)
                                      , clearLegend
                                      , setLegend
                                      , withLegendFormat
                                      -- ** Formatting
                                      , Tick(..), TickValues(..), GridLines
                                      , TickFormat(..)
                                      , setTicks
                                      , setGridlines
                                      , setTickLabelFormat
                                      , setTickLabels
                                      , withTickLabelsFormat
                                      , withAxisLabel
                                      , withAxisLine
                                      , withGridLine
                                       -- * Lines
                                      , Line(), LineFormat()
                                      , DashStyle,Dash(..),LineWidth
                                      , clearLineFormat
                                      , setDashStyle
                                      , setLineWidth
                                      , setLineColour
                                      -- * Points
                                      , Point(), PointFormat()
                                      , Glyph(..)
                                      , PointSize
                                      , setGlyph
                                      , setPointSize
                                      , setPointColour
                                      -- * Bars
                                      , Bar(), BarFormat()
                                      , clearBarFormat
                                      , setBarWidth
                                      , setBarColour
                                      , setBarBorderWidth
                                      , setBarBorderColour
                                      -- * Text labels
                                      , Text()
                                      , FontFamily,FontSize,Color
                                      -- | A text element must exist for formatting to work
                                      , clearText
                                      , clearTextFormat
                                      , setText
                                      , setFontFamily
                                      , setFontStyle
                                      , setFontVariant
                                      , setFontWeight
                                      , setFontStretch
                                      , setFontSize
                                      , setFontColour
                                      ) where

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

--import Data.Packed.Vector
--import Numeric.LinearAlgebra.Linear

--import Data.Word
--import Data.Colour.SRGB
import Data.Colour.Names

import qualified Data.Array.IArray as A

--import qualified Graphics.Rendering.Cairo as C
--import qualified Graphics.Rendering.Pango as P

--import Control.Monad.State
--import Control.Monad.Reader

import Prelude hiding(min,max)

import Graphics.Rendering.Plot.Figure.Text
import Graphics.Rendering.Plot.Figure.Line
import Graphics.Rendering.Plot.Figure.Point
import Graphics.Rendering.Plot.Figure.Bar
import Graphics.Rendering.Plot.Figure.Plot

import Graphics.Rendering.Plot.Types
import Graphics.Rendering.Plot.Defaults

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

-- | perform some actions on the text defaults, must be run before other text element modifications
withTextDefaults :: Text () -> Figure ()
withTextDefaults m = do
  o <- getDefaults
  let to' = _textoptions o
  let (FontText to _) = execText m to' (FontText to' "")
  modifyDefaults $ \s -> s { _textoptions = to }

-- | perform some actions on the line defaults, must be run before other line element modifications
withLineDefaults :: Line () -> Figure ()
withLineDefaults m = do
  o <- getDefaults
  let lo' = _lineoptions o
  let (TypeLine lo _) = execLine m lo' (TypeLine lo' black)
  modifyDefaults $ \s -> s { _lineoptions = lo }

-- | perform some actions on the point defaults, must be run before other point modifications
withPointDefaults :: Point () -> Figure ()
withPointDefaults m = do
  o <- getDefaults
  let po' = _pointoptions o
  let (FullPoint po _) = execPoint m po' (FullPoint po' defaultGlyph)
  modifyDefaults $ \s -> s { _pointoptions = po }

-- | perform some actions on the bar defaults, must be run before other point modifications
withBarDefaults :: Bar () -> Figure ()
withBarDefaults m = do
  o <- getDefaults
  let bo' = _baroptions o
  let (TypeBar bo _) = execBar m bo' (TypeBar bo' black)
  modifyDefaults $ \s -> s { _baroptions = bo }

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

-- | create a new blank 'Figure'
newFigure :: Figure ()
newFigure = putFigure $ Figure defaultFigureBackgroundColour
                          defaultFigurePadding NoText NoText
                          (A.listArray ((1,1),(1,1)) [Nothing])
{-
newLineFigure :: DataSeries                      -- ^ the y series
              -> FigureData
newLineFigure d@(DS_1toN _ _) = let ((xmin,xmax),(ymin,ymax)) = calculateRanges d
                                    plot = Plot False defaultPlotPadding NoText
                                           (defaultRanges xmin xmax ymin ymax)
                                           [defaultXAxis,defaultYAxis]
                                           Nothing Line d []
                                in Figure defaultFigurePadding NoText NoText
                                       (A.listArray ((1,1),(1,1)) [Just plot])
-}
{-
-- | create a new 'Figure'
newFigure :: PlotType -> DataSeries -> Figure ()
newFigure Line   d@(DS_1toN _ _) = putFigure $ newLineFigure d
--newFigure _      _               = error "Figure type not implemented"
-}

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

-- | set the background colour of the figure
setBackgroundColour :: Color -> Figure ()
setBackgroundColour c = modifyFigure $ \s -> s { _back_clr = c }

-- | set the padding of the figure
setFigurePadding :: Double -> Double -> Double -> Double -> Figure ()
setFigurePadding l r b t = modifyFigure $ \s ->
                             s { _fig_pads = Padding l r b t }

-- | operate on the title
withTitle :: Text () -> Figure ()
withTitle m = do
  o <- getDefaults
  modifyFigure $ \s ->
    s { _title = execText m (_textoptions o) (_title s) }

-- | operate on the sub-title
withSubTitle :: Text () -> Figure ()
withSubTitle m = do
  o <- getDefaults
  modifyFigure $ \s ->
    s { _subtitle = execText m (_textoptions o) (_title s) }

-- | set the shape of the plots, losing all current plots
setPlots :: Int      -- ^ rows
         -> Int      -- ^ columns
         -> Figure ()
setPlots r c = modifyFigure $ \s ->
                 s { _plots = A.listArray ((1,1),(r,c))
                      (replicate (r*c) Nothing) }

-- | perform some actions on the specified subplot
withPlot :: (Int,Int) -> Plot () -> Figure ()
withPlot i m = do
  o <- getDefaults
  s <- getSupplies
  modifyFigure $ \p ->
    p { _plots = let plots = _plots p
                     plot' = plots A.! i
                     plot = case plot' of
                              Nothing -> emptyPlot
                              Just p' -> p'
                                                    -- we revert supplies to the original here
                                                    -- since we might want the same colour
                                                    -- order for all plots
                                                    -- HOWEVER: need a better execPlot group
                 in plots A.// [(i,Just $ execPlot m s o plot)] }

-- | perform some actions all subplots
withPlots :: Plot () -> Figure ()
withPlots m = do
  o <- getDefaults
  s <- getSupplies
  modifyFigure $ \p ->
    p { _plots = let plots = _plots p
                     plot p' = case p' of
                                 Nothing  -> emptyPlot
                                 Just p'' -> p''
                 in plots A.// map (\(i,e) ->
                          (i,Just $ execPlot m s o (plot e))) (A.assocs plots) }

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

max_stars_repo_name: JackTheEngineer/plot