text:
{-# OPTIONS_GHC -Wall #-}
{-# LANGUAGE Safe #-}
{-# LANGUAGE MultiParamTypeClasses, FunctionalDependencies #-}
{-# LANGUAGE TypeSynonymInstances, FlexibleInstances #-}
{-# LANGUAGE CPP #-}

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

This module contains ket vectors, bra vectors,
and operators for quantum mechanics.
-}

-- a Ket layer on top of QuantumMat

module Physics.Learn.Ket
    (
    -- * Basic data types
      C
    , i
    , magnitude
    , Ket
    , Bra
    , Operator
    -- * Kets for spin-1/2 particles
    , xp
    , xm
    , yp
    , ym
    , zp
    , zm
    , np
    , nm
    -- * Operators for spin-1/2 particles
    , sx
    , sy
    , sz
    , sn
    , sn'
    -- * Quantum Dynamics
    , timeEvOp
    , timeEv
    -- * Composition
    , Kron(..)
    -- * Measurement
    , possibleOutcomes
    , outcomesProjectors
    , outcomesProbabilities
--    , prob
--    , probs
    -- * Generic multiplication
    , Mult(..)
    -- * Adjoint operation
    , Dagger(..)
    -- * Normalization
    , HasNorm(..)
    -- * Representation
    , Representable(..)
    -- * Orthonormal bases
    , OrthonormalBasis
    , makeOB
    , listBasis
    , size
    -- * Orthonormal bases for spin-1/2 particles
    , xBasis
    , yBasis
    , zBasis
    , nBasis
    -- , angularMomentumXMatrix
    -- , angularMomentumYMatrix
    -- , angularMomentumZMatrix
    -- , angularMomentumPlusMatrix
    -- , angularMomentumMinusMatrix
    -- , jXMatrix
    -- , jYMatrix
    -- , jZMatrix
    -- , matrixCommutator
    -- , rotationMatrix
    -- , jmColumn
    )
    where

-- We try to import only from QuantumMat
-- and not from Numeric.LinearAlgebra

import qualified Data.Complex as C
import Data.Complex
    ( Complex(..)
    , conjugate
    )
import qualified Physics.Learn.QuantumMat as M
import Physics.Learn.QuantumMat
    ( C
    , Vector
    , Matrix
    , (#>)
    , (<#)
    , conjugateTranspose
    , scaleV
    , scaleM
    , conjV
    , fromList
    , toList
    , fromLists
    )
#if MIN_VERSION_base(4,11,0)
import Prelude hiding ((<>))
#endif

infixl 7 <>

-- | A ket vector describes the state of a quantum system.
data Ket = Ket (Vector C)

instance Show Ket where
    show k =
        let message = "Use 'rep <basis name> <ket name>'."
        in if dim k == 2
           then "Representation in zBasis:\n" ++
                show (rep zBasis k) ++ "\n" ++ message
           else message

-- | An operator describes an observable (a Hermitian operator)
--   or an action (a unitary operator).
data Operator = Operator (Matrix C)

instance Show Operator where
    show op =
        let message = "Use 'rep <basis name> <operator name>'."
        in if dim op == 2
           then "Representation in zBasis:\n" ++
                show (rep zBasis op) ++ "\n" ++ message
           else message

-- | A bra vector describes the state of a quantum system.
data Bra = Bra (Vector C)

instance Show Bra where
    show _ = "<bra>\nTry 'rep zBasis <bra name>'"

magnitude :: C -> Double
magnitude = C.magnitude

i :: C
i = 0 :+ 1

-- | Generic multiplication including inner product,
--   outer product, operator product, and whatever else makes sense.
--   No conjugation takes place in this operation.
class Mult a b c | a b -> c where
    (<>) :: a -> b -> c

instance Mult C C C where
    z1 <> z2 = z1 * z2

instance Mult C Ket Ket where
    c <> Ket matrixKet = Ket (scaleV c matrixKet)

instance Mult C Bra Bra where
    c <> Bra matrixBra = Bra (scaleV c matrixBra)

instance Mult C Operator Operator where
    c <> Operator m = Operator (scaleM c m)

instance Mult Ket C Ket where
    Ket matrixKet <> c = Ket (scaleV c matrixKet)

instance Mult Bra C Bra where
    Bra matrixBra <> c = Bra (scaleV c matrixBra)

instance Mult Operator C Operator where
    Operator m <> c = Operator (scaleM c m)

instance Mult Bra Ket C where
    Bra matrixBra <> Ket matrixKet
        = sum $ zipWith (*) (toList matrixBra) (toList matrixKet)

instance Mult Bra Operator Bra where
    Bra matrixBra <> Operator matrixOp
        = Bra (matrixBra <# matrixOp)

instance Mult Operator Ket Ket where
    Operator matrixOp <> Ket matrixKet
        = Ket (matrixOp #> matrixKet)

instance Mult Ket Bra Operator where
    Ket k <> Bra b =
        Operator
        (fromLists [[ x*y | y <- toList b] | x <- toList k])

instance Mult Operator Operator Operator where
    Operator m1 <> Operator m2 = Operator (m1 M.<> m2)

instance Num Ket where
    Ket v1 + Ket v2 = Ket (v1 + v2)
    Ket v1 - Ket v2 = Ket (v1 - v2)
    (*)         = error "Multiplication is not defined for kets"
    negate (Ket v) = Ket (negate v)
    abs         = error "abs is not defined for kets"
    signum      = error "signum is not defined for kets"
    fromInteger = error "fromInteger is not defined for kets"

instance Num Bra where
    Bra v1 + Bra v2 = Bra (v1 + v2)
    Bra v1 - Bra v2 = Bra (v1 - v2)
    (*)         = error "Multiplication is not defined for bra vectors"
    negate (Bra v) = Bra (negate v)
    abs         = error "abs is not defined for bra vectors"
    signum      = error "signum is not defined for bra vectors"
    fromInteger = error "fromInteger is not defined for bra vectors"

instance Num Operator where
    Operator v1 + Operator v2 = Operator (v1 + v2)
    Operator v1 - Operator v2 = Operator (v1 - v2)
    Operator v1 * Operator v2 = Operator (v1 M.<> v2)
    negate (Operator v) = Operator (negate v)
    abs         = error "abs is not defined for operators"
    signum      = error "signum is not defined for operators"
    fromInteger = error "fromInteger is not defined for operators"

-- | The adjoint operation on complex numbers, kets,
--   bras, and operators.
class Dagger a b | a -> b where
    dagger :: a -> b

instance Dagger Ket Bra where
    dagger (Ket v) = Bra (conjV v)

instance Dagger Bra Ket where
    dagger (Bra v) = Ket (conjV v)

instance Dagger Operator Operator where
    dagger (Operator m) = Operator (conjugateTranspose m)

instance Dagger C C where
    dagger c = conjugate c

class HasNorm a where
    norm      :: a -> Double
    normalize :: a -> a

instance HasNorm Ket where
    norm (Ket v) = M.norm v
    normalize k  = (1 / norm k :+ 0) <> k

instance HasNorm Bra where
    norm (Bra v) = M.norm v
    normalize b  = (1 / norm b :+ 0) <> b

-- | An orthonormal basis of kets.
newtype OrthonormalBasis = OB [Ket]
    deriving (Show)

-- | Make an orthonormal basis from a list of linearly independent kets.
makeOB :: [Ket] -> OrthonormalBasis
makeOB = OB . gramSchmidt

size :: OrthonormalBasis -> Int
size (OB ks) = length ks

listBasis :: OrthonormalBasis -> [Ket]
listBasis (OB ks) = ks

{-
newOrthonormalBasis :: Int -> OrthonormalBasis
newOrthonormalBasis = undefined
-}

class Representable a b | a -> b where
    rep :: OrthonormalBasis -> a -> b
    dim :: a -> Int

instance Representable Ket (Vector C) where
    rep (OB ks) k = fromList (map (\bk -> dagger bk <> k) ks)
    dim (Ket v) = M.dim v

instance Representable Bra (Vector C) where
    rep (OB ks) b = fromList (map (\bk -> b <> bk) ks)
    dim (Bra v) = M.dim v

instance Representable Operator (Matrix C) where
    rep (OB ks) op = fromLists [[ dagger k1 <> op <> k2 | k2 <- ks ] | k1 <- ks ]
    dim (Operator m) = let (p,q) = M.size m
                       in if p == q then p else error "dim: non-square operator"

--------------
-- Spin 1/2 --
--------------

-- | State of a spin-1/2 particle if measurement
--   in the x-direction would give angular momentum +hbar/2.
xp :: Ket
xp = Ket M.xp

-- | State of a spin-1/2 particle if measurement
--   in the x-direction would give angular momentum -hbar/2.
xm :: Ket
xm = Ket M.xm

-- | State of a spin-1/2 particle if measurement
--   in the y-direction would give angular momentum +hbar/2.
yp :: Ket
yp = Ket M.yp

-- | State of a spin-1/2 particle if measurement
--   in the y-direction would give angular momentum -hbar/2.
ym :: Ket
ym = Ket M.ym

-- | State of a spin-1/2 particle if measurement
--   in the z-direction would give angular momentum +hbar/2.
zp :: Ket
zp = Ket M.zp

-- | State of a spin-1/2 particle if measurement
--   in the z-direction would give angular momentum -hbar/2.
zm :: Ket
zm = Ket M.zm

-- | State of a spin-1/2 particle if measurement
--   in the n-direction, described by spherical polar angle theta
--   and azimuthal angle phi, would give angular momentum +hbar/2.
np :: Double -> Double -> Ket
np theta phi
    = (cos (theta / 2) :+ 0) <> zp
      + (sin (theta / 2) :+ 0) * (cos phi :+ sin phi) <> zm

-- | State of a spin-1/2 particle if measurement
--   in the n-direction, described by spherical polar angle theta
--   and azimuthal angle phi, would give angular momentum -hbar/2.
nm :: Double -> Double -> Ket
nm theta phi
    = (sin (theta / 2) :+ 0) <> zp
      - (cos (theta / 2) :+ 0) * (cos phi :+ sin phi) <> zm

-- | The orthonormal basis composed of 'xp' and 'xm'.
xBasis :: OrthonormalBasis
xBasis = makeOB [xp,xm]

-- | The orthonormal basis composed of 'yp' and 'ym'.
yBasis :: OrthonormalBasis
yBasis = makeOB [yp,ym]

-- | The orthonormal basis composed of 'zp' and 'zm'.
zBasis :: OrthonormalBasis
zBasis = makeOB [zp,zm]

-- | Given spherical polar angle theta and azimuthal angle phi,
--   the orthonormal basis composed of 'np' theta phi and 'nm' theta phi.
nBasis :: Double -> Double -> OrthonormalBasis
nBasis theta phi = makeOB [np theta phi,nm theta phi]

-- | The Pauli X operator.
sx :: Operator
sx = xp <> dagger xp - xm <> dagger xm

-- | The Pauli Y operator.
sy :: Operator
sy = yp <> dagger yp - ym <> dagger ym

-- | The Pauli Z operator.
sz :: Operator
sz = zp <> dagger zp - zm <> dagger zm

-- | Pauli operator for an arbitrary direction given
--   by spherical coordinates theta and phi.
sn :: Double -> Double -> Operator
sn theta phi
    = (sin theta * cos phi :+ 0) <> sx +
      (sin theta * sin phi :+ 0) <> sy +
      (cos theta           :+ 0) <> sz

-- | Alternative definition
--   of Pauli operator for an arbitrary direction.
sn' :: Double -> Double -> Operator
sn' theta phi
    = np theta phi <> dagger (np theta phi) -
      nm theta phi <> dagger (nm theta phi)

----------------------
-- Quantum Dynamics --
----------------------

-- | Given a time step and a Hamiltonian operator,
--   produce a unitary time evolution operator.
--   Unless you really need the time evolution operator,
--   it is better to use 'timeEv', which gives the
--   same numerical results without doing an explicit
--   matrix inversion.  The function assumes hbar = 1.
timeEvOp :: Double -> Operator -> Operator
timeEvOp dt (Operator m) = Operator (M.timeEvMat dt m)

-- | Given a time step and a Hamiltonian operator,
--   advance the state ket using the Schrodinger equation.
--   This method should be faster than using 'timeEvOp'
--   since it solves a linear system rather than calculating
--   an inverse matrix.  The function assumes hbar = 1.
timeEv :: Double -> Operator -> Ket -> Ket
timeEv dt (Operator m) (Ket k) = Ket $ M.timeEv dt m k

-----------------
-- Composition --
-----------------

class Kron a where
    kron :: a -> a -> a

instance Kron Ket where
    kron (Ket v1) (Ket v2) = Ket (M.kron v1 v2)

instance Kron Bra where
    kron (Bra v1) (Bra v2) = Bra (M.kron v1 v2)

instance Kron Operator where
    kron (Operator m1) (Operator m2) = Operator (M.kron m1 m2)

-----------------
-- Measurement --
-----------------

-- | The possible outcomes of a measurement
--   of an observable.
--   These are the eigenvalues of the operator
--   of the observable.
possibleOutcomes :: Operator -> [Double]
possibleOutcomes (Operator observable) = M.possibleOutcomes observable

-- | Given an obervable, return a list of pairs
--   of possible outcomes and projectors
--   for each outcome.
outcomesProjectors :: Operator -> [(Double,Operator)]
outcomesProjectors (Operator m)
    = [(val,Operator p) | (val,p) <- M.outcomesProjectors m]

-- | Given an observable and a state ket, return a list of pairs
--   of possible outcomes and probabilites
--   for each outcome.
outcomesProbabilities :: Operator -> Ket -> [(Double,Double)]
outcomesProbabilities (Operator m) (Ket v)
    = M.outcomesProbabilities m v

{-
prob :: Ket -> Ket -> Double
prob k1 k2 = magnitude c ** 2
    where
      c = dagger k1 <> k2

probs :: OrthonormalBasis -> Ket -> [Double]
probs (OB ks) k = map (\bk -> let c = dagger bk <> k in magnitude c ** 2) ks
-}

{-
----------------------------------------
-- Angular Momentum of arbitrary size --
----------------------------------------

angularMomentumZMatrix :: Rational -> Matrix Cyclotomic
angularMomentumZMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> if r == c then fromRational (j + 1 - fromIntegral r) else 0)

twoJPlusOne :: Rational -> Maybe Int
twoJPlusOne j
    | j >= 0 && (denominator j == 1 || denominator j == 2)  = Just $ fromIntegral $ numerator (2 * j + 1)
    | otherwise                                             = Nothing

angularMomentumPlusMatrix :: Rational -> Matrix Cyclotomic
angularMomentumPlusMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let mr = j + 1 - fromIntegral r
                                             mc = j + 1 - fromIntegral c
                                         in if mr == mc + 1
                                            then sqrtRat (j*(j+1) - mc*mr)
                                            else 0
                              )

angularMomentumMinusMatrix :: Rational -> Matrix Cyclotomic
angularMomentumMinusMatrix j
    = case twoJPlusOne j of
        Nothing -> error "j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let mr = j + 1 - fromIntegral r
                                             mc = j + 1 - fromIntegral c
                                         in if mr + 1 == mc
                                            then sqrtRat (j*(j+1) - mc*mr)
                                            else 0
                              )

angularMomentumXMatrix :: Rational -> Matrix Cyclotomic
angularMomentumXMatrix j
    = scaleMatrix (1/2) (angularMomentumPlusMatrix j + angularMomentumMinusMatrix j)

angularMomentumYMatrix :: Rational -> Matrix Cyclotomic
angularMomentumYMatrix j
    = scaleMatrix (1/(2*i)) (angularMomentumPlusMatrix j - angularMomentumMinusMatrix j)

jXMatrix :: Rational -> Matrix Cyclotomic
jXMatrix = angularMomentumXMatrix

jYMatrix :: Rational -> Matrix Cyclotomic
jYMatrix = angularMomentumYMatrix

jZMatrix :: Rational -> Matrix Cyclotomic
jZMatrix = angularMomentumZMatrix

matrixCommutator :: Matrix Cyclotomic -> Matrix Cyclotomic -> Matrix Cyclotomic
matrixCommutator m1 m2 = m1 * m2 - m2 * m1

-----------------------
-- Rotation matrices --
-----------------------

r2i :: Rational -> Integer
r2i r
    | denominator r == 1  = numerator r
    | otherwise           = error "r2i:  not an integer"

-- from Sakurai, revised, (3.8.33)
-- beta in degrees
smallDRotElement :: Rational -> Rational -> Rational -> Rational -> Cyclotomic
smallDRotElement j m' m beta
    = sum [parity(k-m+m') * sqrtRat (fact(j+m) * fact(j-m) * fact(j+m') * fact(j-m'))
                   / fromRational (fact(j+m-k) * fact(k) * fact(j-k-m') * fact(k-m+m'))
                         * cosDeg (beta/2) ^ r2i(2*j-2*k+m-m')
                         * sinDeg (beta/2) ^ r2i(2*k-m+m') | k <- [max 0 (m-m') .. min (j+m) (j-m')]]

parity :: Rational -> Cyclotomic
parity = fromIntegral . parityInteger . r2i

-- | (-1)^n, where n might be negative
parityInteger :: Integer -> Integer
parityInteger n
    | odd n      = -1
    | otherwise  =  1

factInteger :: Integer -> Integer
factInteger n
    | n == 0     = 1
    | n >  0     = n * factInteger (n-1)
    | otherwise  = error "factInteger called on negative argument"

fact :: Rational -> Rational
fact = fromIntegral . factInteger . r2i

-- | Rotation matrix elements.
--   From Sakurai, Revised Edition, (3.5.50).
--   The matrix desribes a rotation by gamma about the z axis,
--   followed by a rotation by beta about the y axis,
--   followed by a rotation by alpha about the z axis.
bigDRotElement :: Rational  -- ^ j, a nonnegative integer or half-integer
               -> Rational  -- ^ m', an integer or half-integer index for the row
               -> Rational  -- ^ m, an integer or half-integer index for the column
               -> Rational  -- ^ alpha, in degrees
               -> Rational  -- ^ beta, in degrees
               -> Rational  -- ^ gamma, in degrees
               -> Cyclotomic  -- ^ rotation matrix element
bigDRotElement j m' m alpha beta gamma
    = polarRat 1 (-(m' * alpha + m * gamma) / 360) * smallDRotElement j m' m beta

-- | Rotation matrix for a spin-j particle.
--   The matrix desribes a rotation by gamma about the z axis,
--   followed by a rotation by beta about the y axis,
--   followed by a rotation by alpha about the z axis.
rotationMatrix :: Rational  -- ^ j, a nonnegative integer or half-integer
               -> Rational  -- ^ alpha, in degrees
               -> Rational  -- ^ beta, in degrees
               -> Rational  -- ^ gamma, in degrees
               -> Matrix Cyclotomic  -- ^ rotation matrix
rotationMatrix j alpha beta gamma
    = case twoJPlusOne j of
        Nothing -> error "bigDRotMatrix: j must be a nonnegative integer or half-integer"
        Just d  -> matrix d d (\(r,c) -> let m' = j + 1 - fromIntegral r
                                             m  = j + 1 - fromIntegral c
                                         in bigDRotElement j m' m alpha beta gamma
                              )

----------------------------------
-- Angular Momentum eigenstates --
----------------------------------

jmColumn :: Rational -> Rational -> Matrix Cyclotomic
jmColumn j m
    = case twoJPlusOne j of
        Nothing -> error "bigDRotMatrix: j must be a nonnegative integer or half-integer"
        Just d  -> matrix d 1 (\(r,_) -> let m' = j + 1 - fromIntegral r
                                         in if m == m'
                                            then 1
                                            else 0
                              )
-}

------------------
-- Gram-Schmidt --
------------------

-- | Form an orthonormal list of kets from
--   a list of linearly independent kets.
gramSchmidt :: [Ket] -> [Ket]
gramSchmidt [] = []
gramSchmidt [k] = [normalize k]
gramSchmidt (k:ks) = let nks = gramSchmidt ks
                         nk = normalize (foldl (-) k [w <> dagger w <> k | w <- nks])
                     in nk:nks

-- todo:  Clebsch-Gordon coeffs

max_stars_repo_name: walck/learn-physics

max_issues_repo_name: walck/learn-physics

max_forks_repo_name: walck/learn-physics


text:
module Test where

  import Numeric.LinearAlgebra
  import Data.List
  import Debug.Trace

  import Neural
  import Tool

  {- Network test codes-}

  -- Train a network and run validation and test
  runSuite :: Int -> Double -> [Int] -> ([Example], [Example], [Example]) -> IO ()
  runSuite epoch learningRate hiddenDims (training, validation, test) = do
    initialNetwork <- initNetwork hiddenDims (size $ (fst.head) training)
    let network = train epoch learningRate training initialNetwork
    putStrLn "Validation Phase"
    runTest validation (predict network)
    putStrLn "Test Phase"
    runTest test (predict network)

  -- Run test with given examples
  runTest :: [Example] -> (Vector Double -> Vector Double) -> IO ()
  runTest examples predictFunction =
    do
      putStrLn "Confusion Matrix : "
      disp 6 confusion
      putStr "Precision : "
      putStrLn $ show (getPrecision confusion)
        where
          confusion = getConfusion predictions gts
            where
              predictions = (predictFunction . fst) `map` examples
              gts = snd `map` examples

  -- Generate confusion matrix
  -- Row for prediction, col for ground truth
  getConfusion :: [Vector Double] -> [Vector Double] -> Matrix Double
  --getConfusion pred gts | trace ("getConfusion " ++ show pred ++ " " ++ show gts) False = undefined
  getConfusion pred gts =
    accum (zeromat width height) (+) ((\p -> (p, 1)) `map` pts)
      where
        width = size (head pred)
        height = width
        pts = zip (stepLabel `map` pred) (stepLabel `map` gts)

  -- Computes precision from a confusion matrix
  getPrecision :: Matrix Double -> Double
  getPrecision confusion = (sumElements $ takeDiag confusion) / (sumElements confusion)

  -- Splits examples into test/validation/test sets
  splitDataSet :: [Example] -> ([Example], [Example], [Example])
  splitDataSet wholeData = training `seq` validation `seq` test `seq` (training, validation, test)
    where
      half = length wholeData `quot` 2
      valLen  = 2 * (length wholeData `quot` 10)
      testLen  = 3 * (length wholeData `quot` 10)

      training = take half wholeData
      validation = ((take valLen) . (drop half)) wholeData
      test = ((take testLen) . (drop valLen) . (drop half)) wholeData

ext: hs

max_stars_repo_name: ycjungSubhuman/hskmlp

max_issues_repo_name: ycjungSubhuman/hskmlp

max_forks_repo_name: ycjungSubhuman/hskmlp


text:
{-# LANGUAGE OverloadedLists #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TemplateHaskell #-}
module Myelin.Neuron where

import Control.Lens hiding ((.=), (*~))
import Data.Aeson
import Numeric.LinearAlgebra

{--| NeuronType corresponds to the different point neuron
types that are constructible in SNN. Many of them
are back end specific.

Things that can be improved:

- Currently none of the parameters have any units,
  instead their units are indicated by their names

  - tau_* are time constants
  - v_* are voltages
  - i_* are currents
  - cm is a capacitance

- There is no restriction to "biological" ranges of
  parameters in place and no check of consistency

- In the case of the Heidelberg Hardware system additionally
  every neuron on the chip has analog parameter variation,
  and finite (much coarser than floating point) precision
  in parameter adjustment.

TODO:

HH_cond_exp
GIF_cond_exp

SpikeSourceGamma
SpikeSourceArray
SpikeSourcePoisson
SpikeSourcePoissonRefractory
--}
data NeuronType =
    IFCurrentAlpha {
        _tau_m :: Float, -- ^ membrane time constant
        _tau_refrac :: Float, -- ^ refractory time
        _v_thresh :: Float, -- ^ threshhold potential
        _tau_syn_E :: Float, -- ^ excitatory synaptic time constant
        _v_rest :: Float, -- ^ resting potential
        _cm :: Float, -- ^ membrane capactitance
        _v_reset :: Float, -- ^ reset potential
        _tau_syn_I :: Float, -- ^ inhibitory synaptic time constant
        _i_offset :: Float -- ^ offset current
    }
    | IFCondAlpha {
        _v_rest :: Float, -- ^ resting membrane potential
        _cm :: Float, -- ^ capacity of the membrane
        _tau_m :: Float, -- ^ membrane time constant
        _tau_refrac :: Float, -- ^ duration of refractory period
        _tau_syn_E :: Float, -- ^ rise time of the excitatory synaptic alpha function
        _tau_syn_I :: Float, -- ^ rise time of the inhibitory synaptic alpha function
        _e_rev_E :: Float, -- ^ reversal potential for excitatory input
        _e_rev_I :: Float, -- ^ reversal potential for inhibitory input
        _v_thresh :: Float, -- ^ spike threshold
        _v_reset :: Float, -- ^ reset potential after a spike
        _i_offset :: Float -- ^ offset current
    }
    | IFSpikey {
        _e_rev_I :: Float, -- ^ excitatory reversal current
        _g_leak :: Float, -- ^ leak conductance
        _tau_refrac :: Float, -- ^ refractory time
        _v_reset :: Float, -- ^ reset potential
        _v_rest :: Float, -- ^ resting potential
        _v_thresh :: Float -- ^ threshhold potential
    }
    | IFCurrExp {
        _cm :: Float, -- ^ membrane capacitance
        _tau_m :: Float, -- ^ membrane time constant
        _tau_syn_E :: Float, -- ^ excitatory synaptic time constant
        _tau_syn_I :: Float, -- ^ inhibitory synaptic time constant
        _tau_refrac :: Float, -- ^ refractory time
        _v_thresh :: Float, -- ^ threshhold voltage
        _v_rest :: Float, -- ^ resting potential
        _v_reset :: Float, -- ^ reset potential
        _i_offset :: Float -- ^ offset current
    }
    | IFCondExp {
        _v_rest :: Float, -- ^ resting potential
        _cm :: Float, -- ^ membrane capacitance
        _tau_m :: Float, -- ^ membrane time constant
        _tau_refrac :: Float, -- ^ refractory time
        _tau_syn_E :: Float, -- ^ excitatory synaptic time constant
        _tau_syn_I :: Float, -- ^ inhibitory synaptic time constant
        _e_rev_E :: Float, -- ^ excitatory reversal potential
        _e_rev_I :: Float, -- ^ inhibitory reversal potential
        _v_thresh :: Float, -- ^ spike initiation threshold
        _v_reset :: Float, -- ^ reset value for membrane potential after a spike
        _i_offset :: Float -- ^ offset current
    }
    | Izhikevich {
        _a :: Float,
        _i_offset :: Float,
        _c :: Float,
        _d :: Float,
        _b :: Float
    }
    | EIFCondExp {
        _v_reset :: Float,
        _i_offset :: Float,
        _tau_w :: Float,
        _tau_syn_I :: Float,
        _e_rev_E :: Float,
        _v_rest :: Float,
        _cm :: Float,
        _tau_syn_E :: Float,
        _tau_m :: Float,
        _a :: Float,
        _delta_T :: Float,
        _v_thresh :: Float,
        _b :: Float,
        _v_spike :: Float,
        _e_rev_I :: Float,
        _tau_refrac :: Float
    }
    | EIFCondAlpha {
        _cm :: Float, -- ^ capacity of the membrane
        _tau_refrac :: Float, -- ^ duration of the refractory period
        _v_spike :: Float, -- ^ spike detection threshold
        _v_reset :: Float, -- ^ reset value for membrane potential after a spike
        _v_rest :: Float, -- ^ resting membrane potential (Leak reversal potential)
        _tau_m :: Float, -- ^ membrane time constant
        _i_offset :: Float, -- ^ offset current
        _a :: Float, -- ^ subthreshold adaptation conductance
        _b :: Float, -- ^ spike-triggered adaptation
        _delta_T :: Float, -- ^ slope factor
        _tau_w :: Float, -- ^ adaptation time constant
        _v_thresh :: Float, -- ^ spike initiation threshold
        _e_rev_E :: Float, -- ^ excitatory reversal potential
        _tau_syn_E :: Float, -- ^ rise time of excitatory synaptic conductance (alpha function)
        _e_rev_I :: Float, -- ^ inhibitory reversal potential
        _tau_syn_I :: Float -- ^ rise time of the inhibitory synaptic conductance (alpha function)
    }
    | HHCondExp {
        _cm :: Float, -- ^ capacity of the membrane
        _e_rev_E :: Float, -- ^ excitatory reversal potential
        _e_rev_I :: Float, -- ^ inhibitory reversal potential
        _e_rev_K :: Float,
        _e_rev_Na :: Float,
        _e_rev_leak :: Float,
        _g_leak :: Float,
        _gbar_K :: Float,
        _gbar_Na :: Float,
        _i_offset :: Float,
        _tau_syn_E :: Float,
        _tau_syn_I :: Float,
        _v_offset :: Float
    } deriving (Eq, Show)

makeLenses ''NeuronType
makePrisms ''NeuronType

-- ^  Defaults taken from http://neuralensemble.org/docs/PyNN/standardmodels.html
if_cond_exp :: NeuronType
if_cond_exp = IFCondExp {
    _v_rest = -65.0,
    _cm = 1.0,
    _tau_m = 20.0,
    _tau_refrac = 0.0,
    _tau_syn_E = 5.0,
    _tau_syn_I = 5.0,
    _e_rev_E = 0.0,
    _e_rev_I = -70.0,
    _v_thresh = -50.0,
    _v_reset = -65.0,
    _i_offset = 0.0
}

if_current_alpha :: NeuronType
if_current_alpha = IFCurrentAlpha {
    _tau_m = 20.0,
    _tau_refrac = 0.1,
    _v_thresh = -50.0,
    _tau_syn_E =  0.5,
    _tau_syn_I = 0.5,
    _v_rest = -65.0,
    _cm = 1.0,
    _v_reset = -65.0,
    _i_offset = 0.0
}

if_cond_alpha :: NeuronType
if_cond_alpha = IFCondAlpha {
    _v_rest = -65.0,
    _cm = 1.0,
    _tau_m = 20.0,
    _tau_refrac = 0.0,
    _tau_syn_E = 5.0,
    _tau_syn_I = 5.0,
    _e_rev_E = 0.0,
    _e_rev_I = -70.0,
    _v_thresh = -50.0,
    _v_reset = -65.0,
    _i_offset = 0.0
}

if_spikey :: NeuronType
if_spikey = IFSpikey {
    _e_rev_I = -80.0,
    _g_leak = 20.0,
    _tau_refrac = 1.0,
    _v_reset = -80.0,
    _v_rest = -75.0,
    _v_thresh = -55.0
}

if_current_exponential :: NeuronType
if_current_exponential = IFCurrExp {
    _cm = 1.0,
    _tau_m = 20.0,
    _tau_syn_E = 5.0,
    _tau_syn_I = 5.0,
    _tau_refrac = 0.1,
    _v_thresh = -50.0,
    _v_rest = -65.0,
    _v_reset = -65.0,
    _i_offset = 0.0
}

izhikevich :: NeuronType
izhikevich = Izhikevich {
    _a = 0.02,
    _i_offset = 0.0,
    _c = -65.0,
    _d = 2.0,
    _b = 0.2
}

eif_cond_exp :: NeuronType
eif_cond_exp = EIFCondExp {
    _v_reset = -70.6,
    _i_offset = 0.0,
    _tau_w = 144.0,
    _tau_syn_I = 5.0,
    _e_rev_E = 0.0,
    _v_rest = -70.6,
    _cm = 0.281,
    _tau_syn_E = 5.0,
    _tau_m = 9.3667,
    _a = 4.0,
    _delta_T = 2.0,
    _v_thresh = -50.4,
    _b = 0.0805,
    _v_spike = -40.0,
    _e_rev_I = -80.0,
    _tau_refrac = 0.1
}

eif_cond_alpha :: NeuronType
eif_cond_alpha = EIFCondAlpha {
    _v_reset = -70.6,
    _i_offset = 0.0,
    _tau_w = 144.0,
    _tau_syn_I = 5.0,
    _e_rev_E = 0.0,
    _v_rest = -70.6,
    _cm = 0.281,
    _tau_syn_E = 5.0,
    _tau_m = 9.3667,
    _a = 4.0,
    _delta_T = 2.0,
    _v_thresh = -50.4,
    _b = 0.0805,
    _v_spike = -40.0,
    _e_rev_I = -80.0,
    _tau_refrac = 0.1
}

instance ToJSON NeuronType where
    toJSON IFCondExp {..} = object [
            "type" .= ("IFCondExp" :: String),
            "tau_m" .= _tau_m,
            "tau_refrac" .= _tau_refrac,
            "v_thresh" .= _v_thresh,
            "tau_syn_E" .= _tau_syn_E,
            "v_rest" .= _v_rest,
            "cm" .= _cm,
            "v_reset" .= _v_reset,
            "tau_syn_I" .= _tau_syn_I,
            "i_offset" .= _i_offset,
            "e_rev_E" .= _e_rev_E,
            "e_rev_I" .= _e_rev_I
        ]
    toJSON IFCurrentAlpha {..} = object [
            "type" .= ("IFCurrentAlpha" :: String),
            "tau_m" .= _tau_m,
            "tau_refrac" .= _tau_refrac,
            "v_thresh" .= _v_thresh,
            "tau_syn_E" .= _tau_syn_E,
            "v_rest" .= _v_rest,
            "cm" .= _cm,
            "v_reset" .= _v_reset,
            "tau_syn_I" .= _tau_syn_I,
            "i_offset" .= _i_offset
        ]
    toJSON IFCondAlpha {..} = object [
            "type" .= ("IFCondAlpha" :: String),
            "v_rest" .= _v_rest,
            "cm" .= _cm,
            "tau_m" .= _tau_m,
            "tau_refrac" .= _tau_refrac,
            "tau_syn_E" .= _tau_syn_E,
            "tau_syn_I" .= _tau_syn_I,
            "e_rev_E" .= _e_rev_E,
            "e_rev_I" .= _e_rev_I,
            "v_thresh" .= _v_thresh,
            "v_reset" .= _v_reset,
            "i_offset" .= _i_offset
        ]
    toJSON IFSpikey {..} = object [
            "type" .= ("IFSpikey" :: String),
            "e_rev_I" .= _e_rev_I,
            "g_leak" .= _g_leak,
            "tau_refrac" .= _tau_refrac,
            "v_reset" .= _v_reset,
            "v_rest" .= _v_rest,
            "v_thresh" .= _v_thresh
        ]
    toJSON IFCurrExp {..} = object [
            "type" .= ("IFCurrentExp" :: String),
            "cm" .= _cm,
            "tau_m" .= _tau_m,
            "tau_syn_E" .= _tau_syn_E,
            "tau_syn_I" .= _tau_syn_I,
            "tau_refrac" .= _tau_refrac,
            "v_thresh" .= _v_thresh,
            "v_rest" .= _v_rest,
            "v_reset" .= _v_reset,
            "i_offset" .= _i_offset
        ]
    toJSON Izhikevich {..} = object [
            "type" .= ("Izhikevich" :: String),
            "a" .= _a,
            "i_offset" .= _i_offset,
            "c" .= _c,
            "d" .= _d,
            "b" .= _b
        ]
    toJSON EIFCondExp {..} = object [
            "type" .= ("EIFCondExp" :: String),
            "v_reset" .= _v_reset,
            "i_offset" .= _i_offset,
            "tau_w" .= _tau_w,
            "tau_syn_I" .= _tau_syn_I,
            "e_rev_E" .= _e_rev_E,
            "v_rest" .= _v_rest,
            "cm" .= _cm,
            "tau_syn_E" .= _tau_syn_E,
            "tau_m" .= _tau_m,
            "a" .= _a,
            "delta_T" .= _delta_T,
            "v_thresh" .= _v_thresh,
            "b" .= _b,
            "v_spike" .= _v_spike,
            "e_rev_I" .= _e_rev_I,
            "tau_refrac" .= _tau_refrac
        ]
    toJSON EIFCondAlpha {..} = object [
            "type" .= ("EIFCondAlpha" :: String),
            "cm" .= _cm,
            "tau_refrac" .= _tau_refrac,
            "v_spike" .= _v_spike,
            "v_reset" .= _v_reset,
            "v_rest" .= _v_rest,
            "tau_m" .= _tau_m,
            "i_offset" .= _i_offset,
            "a" .= _a,
            "b" .= _b,
            "delta_T" .= _delta_T,
            "tau_w" .= _tau_w,
            "v_thresh" .= _v_thresh,
            "e_rev_E" .= _e_rev_E,
            "tau_syn_E" .= _tau_syn_E,
            "e_rev_I" .= _e_rev_I,
            "tau_syn_I" .= _tau_syn_I
        ]

instance FromJSON NeuronType where
    parseJSON = withObject "neuron" $ \o -> do
        typ :: String <- o .: "type"
        case typ of
            "IFCondExp" -> IFCondExp <$>
                o .: "v_rest" <*>
                o .: "cm" <*>
                o .: "tau_m" <*>
                o .: "tau_refrac" <*>
                o .: "tau_syn_E" <*>
                o .: "tau_syn_I" <*>
                o .: "e_rev_E" <*>
                o .: "e_rev_I" <*>
                o .: "v_thresh" <*>
                o .: "v_reset" <*>
                o .: "i_offset"
            "IFCurrentAlpha" -> IFCurrentAlpha <$>
                o .: "tau_m" <*>
                o .: "tau_refrac" <*>
                o .: "v_thresh" <*>
                o .: "tau_syn_E" <*>
                o .: "v_rest" <*>
                o .: "cm" <*>
                o .: "v_reset" <*>
                o .: "tau_syn_I" <*>
                o .: "i_offset"
            "IFCondAlpha" -> IFCondAlpha <$>
                o .: "v_rest" <*>
                o .: "cm" <*>
                o .: "tau_m" <*>
                o .: "tau_refrac" <*>
                o .: "tau_syn_E" <*>
                o .: "tau_syn_I" <*>
                o .: "e_rev_E" <*>
                o .: "e_rev_I" <*>
                o .: "v_thresh" <*>
                o .: "v_reset" <*>
                o .: "i_offset"
            "IFSpikey" -> IFSpikey <$>
                o .: "e_rev_I" <*>
                o .: "g_leak" <*>
                o .: "tau_refrac" <*>
                o .: "v_reset" <*>
                o .: "v_rest" <*>
                o .: "v_thresh"
            "IFCurrExp" -> IFCurrExp <$>
                o .: "cm" <*>
                o .: "tau_m" <*>
                o .: "tau_syn_E" <*>
                o .: "tau_syn_I" <*>
                o .: "tau_refrac" <*>
                o .: "v_thresh" <*>
                o .: "v_rest" <*>
                o .: "v_reset" <*>
                o .: "i_offset"
            "Izhikevich" -> Izhikevich <$>
                o .: "a" <*>
                o .: "i_offset"<*>
                o .: "c" <*>
                o .: "d" <*>
                o .: "b"
            "EIFCondExp" -> EIFCondExp <$>
                o .: "v_reset" <*>
                o .: "i_offset" <*>
                o .: "tau_w" <*>
                o .: "tau_syn_I" <*>
                o .: "e_rev_E" <*>
                o .: "v_rest" <*>
                o .: "cm" <*>
                o .: "tau_syn_E" <*>
                o .: "tau_m" <*>
                o .: "a" <*>
                o .: "delta_T" <*>
                o .: "v_thresh" <*>
                o .: "b" <*>
                o .: "v_spike" <*>
                o .: "e_rev_I" <*>
                o .: "tau_refrac"
            "EIFCondAlpha" -> EIFCondAlpha <$>
                o .: "cm" <*>
                o .: "tau_refrac" <*>
                o .: "v_spike" <*>
                o .: "v_reset" <*>
                o .: "v_rest" <*>
                o .: "tau_m" <*>
                o .: "i_offset" <*>
                o .: "a" <*>
                o .: "b" <*>
                o .: "delta_T" <*>
                o .: "tau_w" <*>
                o .: "v_thresh" <*>
                o .: "e_rev_E" <*>
                o .: "tau_syn_E" <*>
                o .: "e_rev_I" <*>
                o .: "tau_syn_I"

max_stars_repo_name: volr/myelin

max_issues_repo_name: volr/myelin

max_forks_repo_name: volr/myelin


text:
{-# LANGUAGE OverloadedLists     #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Isumi.Math.FFT.BruunSpec
  ( spec
  ) where

import           Test.Hspec
import           Test.QuickCheck

import           Data.Complex
import qualified Data.Vector.Unboxed     as UV
import           Isumi.Hspec             (shouldBeAppx)
import           Isumi.Math.AppxEq       ((~=))
import           Isumi.Math.FFT.Bruun
import           Isumi.Math.FFT.SpecUtil

spec :: Spec
spec =
  describe "fftBruun/fftBruunC" $ do
    it "works for a simple sample" $
      fftBruun [1..4] `shouldBeAppx` Just
        [10 :+ 0, (-2) :+ 2, (-2) :+ 0, (-2) :+ (-2)]
    it "works for a simple complex sample" $
      fftBruunC [1 :+ 4, 2 :+ 3, 3 :+ 2, 4 :+ 1] `shouldBeAppx` Just
        [10 :+ 10, 0 :+ 4, (-2) :+ 2, (-4) :+ 0]
    it "works for input with trailing zeros" $
      fftBruun [1, 1, 0, 0] `shouldBeAppx` Just
        [2 :+ 0, 1 :+ (-1), 0, 1 :+ 1]
    it "gives same answer for all inputs as fftDirect" $ property $
      \(VectorOf2ExpSize xs) ->
        Just (fftDirect (realVToComplexV xs)) ~= fftBruun xs
    it "gives same answer for all complex inputs as fftDirect" $ property $
      \(VectorOf2ExpSize rs) (VectorOf2ExpSize cs) ->
        let xs = UV.zipWith (:+) rs cs
         in Just (fftDirect xs) ~= fftBruunC xs

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}

module Section_4_5.TwoLayer where

import Data.List

import Numeric.LinearAlgebra
import System.Random.MWC
import System.Random.MWC.Distributions (normal)

crossEntropyError :: Vector Double -> Vector Double -> Double
crossEntropyError y t =
  let delta = 1.0e-7
   in -1 * (sumElements $ t * cmap (\x -> log (x + delta)) y)

class Indexes c t where
  indexes :: c t -> [IndexOf c]

instance Container Vector t => Indexes Vector t where
  indexes v = [0 .. size v - 1]

instance (Num t, Container Vector t) => Indexes Matrix t where
  indexes m = let (x, y) = size m in (,) <$> [0..x-1] <*> [0..y-1]

numericalGradient :: (Container c Double, Additive (c Double), Indexes c Double) => (c Double -> Double) -> c Double -> c Double
numericalGradient f x =
  foldr add (assoc (size x) 0.0 []) $
    flip map (indexes x) $ \i ->
      let df1 = f (x `add` assoc (size x) 0.0 [(i,  h)])
          df2 = f (x `add` assoc (size x) 0.0 [(i, -h)])
          dfdx = (df1 - df2) / (2*h)
       in assoc (size x) 0.0 [(i, dfdx)]
    where h = 1.0e-4 -- 0.0001

type ActivationFunction = Vector Double -> Vector Double

-- | シグモイド関数
sigmoid :: ActivationFunction
sigmoid = cmap (\x -> 1 / (1 + exp (-x)))

-- | ソフトマックス関数
softmax :: ActivationFunction
softmax a =
  let c = maxElement a
      exp_a = cmap (\x -> exp (x - c)) a
      sum_exp_a = sumElements exp_a
   in cmap (/sum_exp_a) exp_a

type Layer = ( Matrix Double      -- 重み行列
             , Vector Double      -- バイアス
             , ActivationFunction -- 活性化関数
             )
type Network = [Layer]

-- | 学習前のニューラルネットワークを生成する
initNetwork :: [Int]                -- 各層のユニット数
            -> [ActivationFunction] -- 各層の活性化関数
            -> IO Network
initNetwork inputs as =
  withSystemRandom $ \gen ->
    sequence (zipWith (build gen) inouts as)
  where
    inouts :: [(Int, Int)]
    inouts = zip inputs (tail inputs)
    build :: GenIO -> (Int, Int) -> ActivationFunction -> IO Layer
    build gen (x, y) a = do
      w <- (x >< y) <$> (sequence $ replicate (x * y) (normal 0 0.01 gen))
      b <- vector <$> (sequence $ replicate y (normal 0 0.01 gen))
      pure (w, b, a)

-- | 認識（推論）を行う
predict :: Network -> Vector Double -> Vector Double
predict nw x = foldl (\x (w, b, h) -> h (x <# w + b)) x nw

-- | 損失関数の値を求める
loss :: [Vector Double] -> [Vector Double] -> Network -> Double
loss xs ts nw =
  let batchSize = genericLength xs
      ys = map (predict nw) xs
   in sum (zipWith crossEntropyError ys ts) / batchSize

-- | 認識精度を求める
accuracy :: [Vector Double] -> [Vector Double] -> Network -> Double
accuracy xs ts nw =
  let batchSize = genericLength xs
      ys = map (predict nw) xs
      yis = map maxIndex ys
      tis = map maxIndex ts
      correctCount = genericLength (filter id (zipWith (==) yis tis))
   in correctCount / batchSize

dismantle :: [a] -> [(a, a -> [a])]
dismantle xs = map (\i -> (xs!!i, (\x -> (take i xs) ++ [x] ++ (drop (i+1) xs)))) [0..length xs - 1] -- TODO more efficient implementation

updateNetwork :: ([Vector Double] -> [Vector Double] -> Network -> Double) -> [Vector Double] -> [Vector Double] -> Network -> Network
updateNetwork loss xs ts network = map (uncurry $ updateLayer xs ts) $ dismantle network
  where
    rate :: Double
    rate = 0.1
    updateLayer :: [Vector Double] -> [Vector Double] -> Layer -> (Layer -> Network) -> Layer
    updateLayer xs ts (w, b, a) mkNw =
      let dw = numericalGradient (\w -> loss xs ts (mkNw (w, b, a))) w
          db = numericalGradient (\b -> loss xs ts (mkNw (w, b, a))) b
       in (w - rate `scale` dw, b - rate `scale` db, a)


text:
module Numbers where

import Control.Monad.Except (throwError)
import Data.Complex
import Data.Ratio
import DataTypes
       (Arity(..), LispError(..), LispVal(..), PrimitiveFunc, ThrowsError)

numPrimitives :: [(String, PrimitiveFunc)]
numPrimitives =
  [ ("+", numAdd)
  , ("-", numSub)
  , ("*", numMul)
  , ("/", numDiv)
  , ("modulo", numMod)
  , ("Integer?", isNumber)
  , ("complex?", isComplex)
  , ("real?", isReal)
  , ("rational?", isRational)
  , ("integer?", isInteger)
  , ("=", numBoolBinopEq)
  , ("/=", numBoolBinopNeq)
  , (">", numBoolBinopGt)
  , ("<", numBoolBinopLt)
  , (">=", numBoolBinopGte)
  , ("<=", numBoolBinopLte)
  , ("quotient", numQuotient)
  , ("remainder", numRem)
  , ("sin", numSine)
  , ("cos", numCos)
  ]

foldlM :: Monad m => (a -> b -> m a) -> a -> [b] -> m a
foldlM f v (x:xs) = (f v x) >>= \a -> foldlM f a xs
foldlM _ v [] = return v

foldl1M :: Monad m => (a -> a -> m a) -> [a] -> m a
foldl1M f (x:xs) = foldlM f x xs
foldl1M _ _ = error "Unexpected error in foldl1M"

numAdd :: PrimitiveFunc
numAdd [] = return $ Integer 0
numAdd a = foldlM (\b c -> doAdd =<< (numCast [b, c])) (Integer 0) a
  where
    doAdd :: LispVal -> ThrowsError LispVal
    doAdd (List [Integer c, Integer d]) = return $ Integer (c + d)
    doAdd (List [Rational c, Rational d]) = return $ Rational (c + d)
    doAdd (List [Float c, Float d]) = return $ Float (c + d)
    doAdd (List [Complex c, Complex d]) = return $ Complex (c + d)
    doAdd _ = throwError $ Default "Unexpected error in +"

numSub :: PrimitiveFunc
numSub [] = throwError $ NumArgs (Min 1) 0 []
numSub a = foldl1M (\b c -> doSub =<< (numCast [b, c])) a -- TODO zero args check
  where
    doSub :: LispVal -> ThrowsError LispVal
    doSub (List [Integer c, Integer d]) = return $ Integer (c - d)
    doSub (List [Rational c, Rational d]) = return $ Rational (c - d)
    doSub (List [Float c, Float d]) = return $ Float (c - d)
    doSub (List [Complex c, Complex d]) = return $ Complex (c - d)
    doSub _ = throwError $ Default "Unexpected error in -"

numMul :: PrimitiveFunc
numMul [] = return $ Integer 1
numMul a = foldl1M (\b c -> doMul =<< (numCast [b, c])) a
  where
    doMul :: LispVal -> ThrowsError LispVal
    doMul (List [Integer c, Integer d]) = return $ Integer (c * d)
    doMul (List [Rational c, Rational d]) = return $ Rational (c * d)
    doMul (List [Float c, Float d]) = return $ Float (c * d)
    doMul (List [Complex c, Complex d]) = return $ Complex (c * d)
    doMul _ = throwError $ Default "Unexpected error in -"

numDiv :: PrimitiveFunc
numDiv [] = throwError $ NumArgs (Min 1) 0 []
numDiv a = foldl1M (\b c -> doDiv =<< (numCast [b, c])) a -- TODO: Zero division error
  where
    doDiv :: LispVal -> ThrowsError LispVal
    doDiv (List [Integer c, Integer d]) =
      if d == 0
        then throwError $ Default "Zero division error"
        else return $ Rational (c % d)
    doDiv (List [Rational c, Rational d]) =
      if d == 0
        then throwError $ Default "Zero division error"
        else return $ Rational (c / d)
    doDiv (List [Float c, Float d]) =
      if d == 0.0
        then throwError $ Default "Zero division error"
        else return $ Float (c / d)
    doDiv (List [Complex c, Complex d]) =
      if d == 0
        then throwError $ Default "Zero division error"
        else return $ Complex (c / d)
    doDiv _ = throwError $ Default "Unexpected error in /"

numMod :: PrimitiveFunc
numMod [] = throwError $ NumArgs (MinMax 2 2) 0 []
numMod (a:[]) = throwError $ NumArgs (MinMax 2 2) 1 [a]
numMod [a, b] = do
  c <- numCast $ [a, b]
  doMod c
  where
    doMod :: LispVal -> ThrowsError LispVal
    doMod (List [Integer c, Integer d]) = return $ Integer (c `mod` d)
    doMod (List [c@(Rational _), d@(Rational _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Rational ((c' `mod` d') % 1)
    doMod (List [c@(Float _), d@(Float _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Float (fromInteger (c' `mod` d'))
    doMod (List [c@(Complex _), d@(Complex _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Complex (fromInteger (c' `mod` d') :+ 0)
    doMod _ = throwError $ Default "Unexpected error in modulo"
numMod a = throwError $ NumArgs (MinMax 2 2) (length a) a

numRem :: PrimitiveFunc
numRem [] = throwError $ NumArgs (MinMax 2 2) 0 []
numRem (a:[]) = throwError $ NumArgs (MinMax 2 2) 1 [a]
numRem [a, b] = do
  c <- numCast $ [a, b]
  doRem c
  where
    doRem :: LispVal -> ThrowsError LispVal
    doRem (List [Integer c, Integer d]) = return $ Integer (c `rem` d)
    doRem (List [c@(Rational _), d@(Rational _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Rational ((c' `rem` d') % 1)
    doRem (List [c@(Float _), d@(Float _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Float (fromInteger (c' `rem` d'))
    doRem (List [c@(Complex _), d@(Complex _)]) = do
      Integer c' <- numToInt c
      Integer d' <- numToInt d
      return $ Complex (fromInteger (c' `rem` d') :+ 0)
    doRem _ = throwError $ Default "Unexpected error in remainder"
numRem a = throwError $ NumArgs (MinMax 2 2) (length a) a

isInteger :: PrimitiveFunc
isInteger [Integer _] = return $ Bool True
isInteger [_] = return $ Bool False
isInteger a = throwError $ NumArgs (MinMax 1 1) (length a) a

isRational :: PrimitiveFunc
isRational [Rational _] = return $ Bool True
isRational a = isInteger a

isReal :: PrimitiveFunc
isReal [Float _] = return $ Bool True
isReal a = isRational a

isComplex :: PrimitiveFunc
isComplex [Complex _] = return $ Bool True
isComplex a = isReal a

isNumber :: PrimitiveFunc
isNumber = isComplex

numToInt :: LispVal -> ThrowsError LispVal
numToInt a@(Integer _) = return a
numToInt (Rational a) =
  if denominator a == 1
    then return $ Integer (numerator a)
    else throwError $ Default "Could not convert rational to integer"
numToInt (Float a) =
  if a == fromInteger (round a)
    then return $ Integer (round a)
    else throwError $ Default "Could not convert float to integer"
numToInt (Complex a) =
  let rp = realPart a
  in if imagPart a == 0 && rp == fromInteger (round rp)
       then return $ Integer (round rp)
       else throwError $ Default "Could not convert complex to integer"
numToInt _ = throwError $ Default "Could not convert non-number to integer"

numCast :: [LispVal] -> ThrowsError LispVal
numCast [a@(Integer _), b@(Integer _)] = return $ List [a, b]
numCast [a@(Rational _), b@(Rational _)] = return $ List [a, b]
numCast [a@(Float _), b@(Float _)] = return $ List [a, b]
numCast [a@(Complex _), b@(Complex _)] = return $ List [a, b]
-- Integer
numCast [(Integer a), b@(Rational _)] = return $ List [Rational (a % 1), b]
numCast [(Integer a), b@(Float _)] = return $ List [Float (fromInteger a), b]
numCast [(Integer a), b@(Complex _)] =
  return $ List [Complex (fromInteger a :+ 0), b]
-- Rational
numCast [a@(Rational _), (Integer b)] = return $ List [a, Rational (b % 1)]
numCast [(Rational a), b@(Float _)] = return $ List [Float (fromRational a), b]
numCast [(Rational a), b@(Complex _)] =
  return $ List [Complex (fromRational a :+ 0), b]
-- Float
numCast [a@(Float _), (Rational b)] = return $ List [a, Float (fromRational b)]
numCast [a@(Float _), (Integer b)] = return $ List [a, Float (fromInteger b)]
numCast [(Float a), b@(Complex _)] = return $ List [Complex (a :+ 0), b]
-- Complex
numCast [a@(Complex _), (Rational b)] =
  return $ List [a, Complex (fromRational b :+ 0)]
numCast [a@(Complex _), (Float b)] = return $ List [a, Complex (b :+ 0)]
numCast [a@(Complex _), (Integer b)] =
  return $ List [a, Complex (fromInteger b :+ 0)]
numCast [a, b] =
  case a of
    Integer _ -> doThrowError b
    Float _ -> doThrowError b
    Rational _ -> doThrowError b
    Complex _ -> doThrowError b
    _ -> doThrowError a
  where
    doThrowError :: LispVal -> ThrowsError LispVal
    doThrowError num = throwError $ TypeMismatch "Integer" num
numCast _ = throwError $ Default "Unexpected error in numCast"

numBoolBinop ::
     String
  -> (LispVal -> LispVal -> ThrowsError LispVal)
  -> LispVal
  -> [LispVal]
  -> ThrowsError LispVal
numBoolBinop name' op b (c:d) = do
  List [b', c'] <- numCast [b, c]
  result <- op b' c'
  case result of
    Bool True -> numBoolBinop name' op c' d
    Bool False -> return $ Bool False
    _ -> throwError $ Default $ "Unexpected error in " ++ name'
numBoolBinop _ _ _ _ = return $ Bool True

numBoolBinopEq :: [LispVal] -> ThrowsError LispVal
numBoolBinopEq [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopEq (x:xs) = numBoolBinop "=" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c == d)
    fn (Rational c) (Rational d) = return $ Bool (c == d)
    fn (Float c) (Float d) = return $ Bool (c == d)
    fn (Complex c) (Complex d) = return $ Bool (c == d)
    fn _ _ = throwError $ Default "Unexpected error in ="

numBoolBinopNeq :: [LispVal] -> ThrowsError LispVal
numBoolBinopNeq [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopNeq (x:xs) = numBoolBinop "/=" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c /= d)
    fn (Rational c) (Rational d) = return $ Bool (c /= d)
    fn (Float c) (Float d) = return $ Bool (c /= d)
    fn (Complex c) (Complex d) = return $ Bool (c /= d)
    fn _ _ = throwError $ Default "Unexpected error in /="

numBoolBinopLt :: [LispVal] -> ThrowsError LispVal
numBoolBinopLt [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopLt (x:xs) = numBoolBinop "<" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c < d)
    fn (Rational c) (Rational d) = return $ Bool (c < d)
    fn (Float c) (Float d) = return $ Bool (c < d)
    fn (Complex _) (Complex _) =
      throwError $ Default "< not defined for complex numbers"
    fn _ _ = throwError $ Default "Unexpected error in <"

numBoolBinopLte :: [LispVal] -> ThrowsError LispVal
numBoolBinopLte [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopLte (x:xs) = numBoolBinop "<=" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c <= d)
    fn (Rational c) (Rational d) = return $ Bool (c <= d)
    fn (Float c) (Float d) = return $ Bool (c <= d)
    fn (Complex _) (Complex _) =
      throwError $ Default "<= not defined for complex numbers"
    fn _ _ = throwError $ Default "Unexpected error in <="

numBoolBinopGt :: [LispVal] -> ThrowsError LispVal
numBoolBinopGt [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopGt (x:xs) = numBoolBinop ">" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c > d)
    fn (Rational c) (Rational d) = return $ Bool (c > d)
    fn (Float c) (Float d) = return $ Bool (c > d)
    fn (Complex _) (Complex _) =
      throwError $ Default "> not defined for complex numbers"
    fn _ _ = throwError $ Default "Unexpected error in >"

numBoolBinopGte :: [LispVal] -> ThrowsError LispVal
numBoolBinopGte [] = throwError $ NumArgs (Min 1) 0 []
numBoolBinopGte (x:xs) = numBoolBinop ">=" fn x xs
  where
    fn :: LispVal -> LispVal -> ThrowsError LispVal
    fn (Integer c) (Integer d) = return $ Bool (c >= d)
    fn (Rational c) (Rational d) = return $ Bool (c >= d)
    fn (Float c) (Float d) = return $ Bool (c >= d)
    fn (Complex _) (Complex _) =
      throwError $ Default ">= not defined for complex numbers"
    fn _ _ = throwError $ Default "Unexpected error in >="

numQuotient :: PrimitiveFunc
numQuotient args =
  if length args /= 2
    then throwError $ NumArgs (MinMax 2 2) (length args) args
    else do
      nums <- numCast args
      case nums of
        List [Integer a, Integer b] -> return $ Integer (a `quot` b)
        _ -> throwError $ Default "Unexpected error in <=" -- TODO better errors

unaryTrig ::
     (Double -> Double) -> (Double -> Double -> Complex Double) -> PrimitiveFunc
unaryTrig op complexOp args =
  if length args /= 1
    then throwError $ NumArgs (MinMax 1 1) (length args) args
    else case args of
           [Integer a] -> return $ Float (op $ fromInteger a)
           [Rational a] -> return $ Float (op $ fromRational a)
           [Float a] -> return $ Float (op a)
           [Complex a] -> return $ Complex $ complexOp (realPart a) (imagPart a)
           _ -> throwError $ Default "Numerical input expected"

numSine :: PrimitiveFunc
numSine = unaryTrig sin (\r i -> ((sin r) * (cosh i)) :+ ((cos r) * (sinh i)))

numCos :: PrimitiveFunc
numCos =
  unaryTrig cos (\r i -> ((cos r) * (cosh i)) :+ (-1 * ((sin r) * (sinh i))))

max_stars_repo_name: ebenpack/hascheme

max_issues_repo_name: ebenpack/hascheme

max_forks_repo_name: ebenpack/hascheme


text:
{-# LANGUAGE BangPatterns #-}
module LINreg2 (mmlLinearRegLP, diagl, rowM, colM,(<+>), (<->), foldl'Rnf
                , mkZ,realignedNpieces,lineFun, mmlCorr-- ,vmean,vvar,vstDiv,vvarM,vunStdNorm,vstNorm, corrRoundCentre
                , standardizedFit,relFit, corr,nCorrV,logg
                ) where

import Numeric.LinearAlgebra
import Numeric.GSL.Statistics as S
import Control.DeepSeq (NFData, rnf)
import Numeric.GSL.Statistics (mean,variance,variance_m,stddev)
import Control.Parallel (pseq)
import Data.Ord (comparing)
import Data.List (transpose, nub,minimumBy,foldl')
import LMisc
import Data.Time
import IO

ident n = (diag . fromList) (take n $ repeat 1) :: Matrix Double

vector xs = fromList xs :: Vector Double
diagl = diag . vector
rowM = asRow . vector
colM = asColumn . vector
-------------------------------------------------------
foldl'Rnf :: (NFData a) => (a -> b -> a) -> a -> [b] -> a
foldl'Rnf f z xs = lgo z xs
    where
        lgo z   []      = z
        lgo z (x:xs)    = lgo z' xs
            where
                z' = rnf (f z x) `pseq` (f z x)

------ standardise by subtracting the mean and dividing by the standard deviation
vstNorm       xs            = mapVector (\ a -> (a - m) / std) xs
                 where m    = mean xs
                       std  = sqrt $ variance_m m xs

-- unstandardising
vunStdNorm m std vs  = mapVector (\a -> std * a + m)  vs

--mean,variance,variance_m,stddev)
vunStdNorm' vs = vunStdNorm (mean vs) (stddev vs)  vs
-----------------------------------------------------------------------------------------------
logg a      --     = log a
      | a <= 0     = a
      | otherwise  = log a
-------------------------------------------------------
infix 1  <+>
(<+>) :: Matrix Double -> Matrix Double -> Matrix Double
m1  <+>    m2   = fromRows $ zipWith (zipVectorWith (+)) (toRows m1) (toRows m2)

(<->) :: Matrix Double -> Matrix Double -> Matrix Double
m1  <->    m2  =  fromRows $ zipWith (zipVectorWith (-)) (toRows m1) (toRows m2)
--}

-- adds a to the diadonal elements of ma
mkZ :: Double -> Matrix Double -> Matrix Double
mkZ  a   ma  = mapMatrixWithIndex (pli a) (trans ma <>  ma)
    where  pli b (!i,!j) v = if i == j then v + b else v

mmlLinearRegLP :: Vector Double ->  Double -> Double -> Matrix Double -> (Double, (Vector Double, (Double, Vector Double))) -- Vector Double
mmlLinearRegLP     xs     m   r  ys  =  (ans, (xTy, (sig, abar)))
                where
                  z       =   mkZ m ys -- make the Z matrix
                  n       =   fromIntegral $ rows  ys
                  !mM     =   n + k  + 1
                  !rest   =   log r + 0.5 * log (2 * mM * abs (det z)) + (n / 2) -- cd (k+1)
                  --abar   =   (inv z) <> (trans ys) <>  xs --
                  abar_aux    =    toColumns $ linearSolveLS z (asColumn ((trans ys) <>  xs))
                  abar    =  head abar_aux
                  !k      =   fromIntegral $ cols z
                  sigSq   =    (xs <.> xs)  - ( xs <.> (ys <> abar))
                  !sig    =   sqrt (sigSq / n)
                  xTy     =     (ys <>  abar) -- predictions fromList . manyInvNS1 Nothing $ toList
                  ans     =   (mM/2) * log (2 * pi) + n * log sig - 0.5 * (k + 1) * log k + rest -- k is 2 here
----
mmlMLwithAbar ::Vector Double ->  Vector Double ->  Double -> Double -> Matrix Double -> Double
mmlMLwithAbar abar     xs     m   r  ys  =  ans
                where
                  z       =   mkZ m ys -- make the Z matrix
                  n       =   fromIntegral $ rows  ys
                  !mM     =   n + k  + 1
                  !k      =   fromIntegral $ cols z
                  !rest   =   log r + 0.5 * logg (2 * mM * abs (det z)) + (n / 2) -- cd (k+1)
                  sigSq   =   (xs <.> xs)  - ( xs <.> (ys <> abar))
                  !sig    =   sqrt (sigSq / n)
                  xTy     =     (ys <>  abar) -- predictions fromList . manyInvNS1 Nothing $ toList
                  ans     =   (mM/2) * log (2 * pi) + n * log sig - 0.5 * (k + 1) * log m + rest

--- realigning a matrix
--- map mkMatrix : [[(Matrix Double,[Int])]]
realignedNpieces :: Int ->  --- number of interventions
    Int ->                  --- minimum length of each sub group
   -- Int ->                  -- maximum length of each subgroup
    [[Double]] ->             --- the list of values
    [(Matrix Double,[Int])] --- the matrix for analysis, along with the locations of interventions
realignedNpieces n m  =  nub . concatMap combMatrix . transpose . map  mkMatrix -- map (\ks -> (newMX, ks)) ctts
    where
        ---
        combMatrix :: [(Matrix Double,[Int])] -> [(Matrix Double, [Int]) ]
        combMatrix mcI  = map (\ks -> (newMX, ks)) cts
            where
                (mxx, cts) = unzip mcI
                newMX      = appendByColumns mxx
        mkMatrix :: [Double] ->  [(Matrix Double,[Int])]
        mkMatrix xs      =   [(reAlign ps ,ks) |
                                let nPnms  =  nPiecesP n m Nothing xs
                                , (ps,ks) <-   nPnms
                                ]
            where
                nPns             = []
                len              = length xs
                cons a vs        = join [fromList [a] , vs]
                app2rows f       = fromRows . map f  . toRows
                tS'              = app2rows (cons 1) . trans . app2rows vstNorm . fromLists --vstNorm
                reAlign          = tS' . adjustMatrix . fillZeroes_aux len 0 1
appendByColumns :: Element t => [Matrix t] -> Matrix t
appendByColumns  = fromRows . concatMap toRows

-- draw a best fit line through a list of points
lineFun :: Vector Double -> (Double,Double)
lineFun vs = (ks @> 1, ks @> 0)
    where
        mx          =  (fst . head . realignedNpieces 0 1 . (: [])) [1 .. fromIntegral $ dim  vs]
        (_,(_,ks))  =  snd $ mmlLinearRegLP (vstNorm vs) 1 1  mx

--
-- finding mormalized correlation using  using MML
mmlCorr :: [Double] -> [Double] -> Double
mmlCorr xs ys =  nCorrV (fromList xs) (fromList ys)

mmlCorr1  xs ys =  (mmlCorr_aux (Just n) xs ys)/ (mmlCorr_aux Nothing [1 .. fromIntegral n] [1 .. fromIntegral n])
    where
        n           = (minimum . map length) [xs, ys]
        mmlCorr_aux mN xs' ys' = (ks @> 1)
            where
                stNorm      =  toList . vstNorm . fromList
                yss         =  (stNorm . maybe ys' (flip take ys')) mN
                vss         =  (vstNorm . fromList . maybe ys' (flip take xs')) mN
                (_,(_,ks))  =  snd $ mmlLinearRegLP vss 1 1  mx
                mx          =  (fst . head . realignedNpieces 0 1 . (: [])) yss

-- relative fit of twoe straight line functions
relFit :: [Double] -> [Double] -> Double
relFit xs ys   = relFit_1 (vstNorm $ fromList xs) (vstNorm $ fromList ys)
relFit_1  xs ys =   abs ((mly  - mlx)  /  mlx)
    where
        mly    = mmlMLwithAbar ab  ys  1 1  (mx ys)
        (mlx,(_,(_,ab))) =   mmlFun  xs
        mx   xs'     =  (fst . head . realignedNpieces 0 1 . (: [])) [1 .. fromIntegral $ dim xs']
        mmlFun vs    =  mmlLinearRegLP vs 1 1  (mx vs)

---pearsons correlation
corr xs  ys  =  S.correlation (fromList xs) (fromList ys)

-- normalised correlation
nCorrV xs ys  = (S.correlation (subVector 0 n xs) (subVector 0 n ys))
    where
         n    =  min (dim xs) (dim ys)

---
standardizedFit  :: [Vector Double] -> Double
standardizedFit [] = 0
standardizedFit ks@(_:_) = S.stddev . fromList $ map mSgLn ks
    where
        mSgLn  ys    = mmlMLwithAbar ab  ys  1 1  (mx ys)
        (_,(_,(_,ab)))   = minimumBy (comparing fst) $ map mmlFun ks -- (fst . snd . snd)
        mx   xs'     =  (fst . head . realignedNpieces 0 1 . (: [])) [1 .. fromIntegral $ dim xs']
        mmlFun vs    =  mmlLinearRegLP vs 1 1  (mx vs)

max_stars_repo_name: rawlep/MML

max_issues_repo_name: rawlep/MML

max_forks_repo_name: rawlep/MML


text:
import Debug.Trace
-- import Control.Monad.State.Lazy

import System.Environment (getArgs)
import System.IO (readFile)
import qualified Data.Set as S
import qualified Data.Map.Strict as M
import Data.Set (Set)
import Data.Map.Strict (Map)
import Data.Tuple (swap)
import Data.Complex (Complex((:+)), polar)

import Data.List (sort, nub)

readMap :: [[Char]] -> [(Int, Int)]
readMap rows = [(x, y) | (y, row) <- zip [0 ..] rows
                       , (x, c) <- zip [0 ..] row
                       , c == '#']

analyzeMap :: [(Int, Int)] -> (Int, (Int, Int))
analyzeMap points = maximum [(length $ nub angles, p)
                            | p <- points
                            , let angles = map (fst . polarDirection p) points]

polarDirection :: (Int, Int) -> (Int, Int) -> (Float, Float)
polarDirection origin@(x0,y0) pt@(x1,y1) = -- (angle_from_vertical, distance)
  swap -- angle first
  $ positiveAngle -- polar's angle is from -pi to pi, convert to 0 to 2pi so that we can sort by angle
  $ polar -- convert to (magnitude, phase)
  $ fromIntegral (-1 * y) :+ fromIntegral x -- make complex, normally
           -- the positive x axis is phase 0, but here we want to go
           -- clockwise from vertical, so swap x and y and change y's
           -- sign to correct for the fact that the coordinates
           -- increase in the down direction.
  where (x, y) = (x1-x0,y1-y0) -- point from perspective of origin
        positiveAngle (m,p) | p < 0 = (m, 2 * pi + p)
                            | otherwise = (m, p)

main = do
  [mapFile] <- getArgs
  m <- readFile mapFile
  let points = readMap $ lines m
  let (reachable, station) = analyzeMap points
  putStrLn $ "Part 1: " ++ show (reachable, station)
  let polarPoints = -- traceShowId $
        sort -- [(round, angle, xy)] in order
        $ foldl (\l ((angle, distance), xy) ->
                   (if null l then 0
                    else let (round0, angle0, pt0) = head l
                         in if angle == angle0 then round0 + 1
                            else 0,
                     angle, xy):l) ([]::[(Int, Float, Int)])
        -- $ traceShowId
        $ sort -- [((angle, distance), xy)] in order
        $ map (\(x,y) -> (polarDirection station (x,y),100*x+y))
        $ filter (/= station) points

  putStrLn $ "Part 2: " ++ show (polarPoints!!199)
  --putStrLn $ "Part 2: " ++ (show $ head $ drop 199 polarPoints)


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE ScopedTypeVariables #-}
#if __GLASGOW_HASKELL__ >= 800
  {-# OPTIONS_GHC -Wno-redundant-constraints #-}
#endif
-- |
-- Module      : Graphics.ColorSpace.Elevator
-- Copyright   : (c) Alexey Kuleshevich 2018-2019
-- License     : BSD3
-- Maintainer  : Alexey Kuleshevich <lehins@yandex.ru>
-- Stability   : experimental
-- Portability : non-portable
--
module Graphics.ColorSpace.Elevator (
  Elevator(..)
  , clamp01
  ) where

import qualified Data.Complex as C
import Data.Int
import Data.Typeable
import Data.Vector.Storable (Storable)
import Data.Vector.Unboxed (Unbox)
import Data.Word
import GHC.Float

-- | A class with a set of convenient functions that allow for changing precision of
-- channels within pixels, while scaling the values to keep them in an appropriate range.
--
-- >>> let rgb = PixelRGB 0.0 0.5 1.0 :: Pixel RGB Double
-- >>> eToWord8 <$> rgb
-- <RGB:(0|128|255)>
-- >>> eToWord16 <$> rgb
-- <RGB:(0|32768|65535)>
--
class (Eq e, Num e, Typeable e, Unbox e, Storable e) => Elevator e where

  -- | Values are scaled to @[0, 255]@ range.
  eToWord8 :: e -> Word8

  -- | Values are scaled to @[0, 65535]@ range.
  eToWord16 :: e -> Word16

  -- | Values are scaled to @[0, 4294967295]@ range.
  eToWord32 :: e -> Word32

  -- | Values are scaled to @[0, 18446744073709551615]@ range.
  eToWord64 :: e -> Word64

  -- | Values are scaled to @[0.0, 1.0]@ range.
  eToFloat :: e -> Float

  -- | Values are scaled to @[0.0, 1.0]@ range.
  eToDouble :: e -> Double

  -- | Values are scaled from @[0.0, 1.0]@ range.
  eFromDouble :: Double -> e

-- | Lower the precision
dropDown :: forall a b. (Integral a, Bounded a, Integral b, Bounded b) => a -> b
dropDown !e = fromIntegral $ fromIntegral e `div` ((maxBound :: a) `div`
                                                   fromIntegral (maxBound :: b))
{-# INLINE dropDown #-}

-- | Increase the precision
raiseUp :: forall a b. (Integral a, Bounded a, Integral b, Bounded b) => a -> b
raiseUp !e = fromIntegral e * ((maxBound :: b) `div` fromIntegral (maxBound :: a))
{-# INLINE raiseUp #-}

-- | Convert to fractional with value less than or equal to 1.
squashTo1 :: forall a b. (Fractional b, Integral a, Bounded a) => a -> b
squashTo1 !e = fromIntegral e / fromIntegral (maxBound :: a)
{-# INLINE squashTo1 #-}

-- | Convert to integral streaching it's value up to a maximum value.
stretch :: forall a b. (RealFrac a, Floating a, Integral b, Bounded b) => a -> b
stretch !e = round (fromIntegral (maxBound :: b) * clamp01 e)
{-# INLINE stretch #-}

-- | Clamp a value to @[0, 1]@ range.
clamp01 :: (Ord a, Floating a) => a -> a
clamp01 !x = min (max 0 x) 1
{-# INLINE clamp01 #-}

-- | Values between @[0, 255]]@
instance Elevator Word8 where
  eToWord8 = id
  {-# INLINE eToWord8 #-}
  eToWord16 = raiseUp
  {-# INLINE eToWord16 #-}
  eToWord32 = raiseUp
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1
  {-# INLINE eToDouble #-}
  eFromDouble = eToWord8
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 65535]]@
instance Elevator Word16 where
  eToWord8 = dropDown
  {-# INLINE eToWord8 #-}
  eToWord16 = id
  {-# INLINE eToWord16 #-}
  eToWord32 = raiseUp
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1
  {-# INLINE eToDouble #-}
  eFromDouble = eToWord16
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 4294967295]@
instance Elevator Word32 where
  eToWord8 = dropDown
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown
  {-# INLINE eToWord16 #-}
  eToWord32 = id
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1
  {-# INLINE eToDouble #-}
  eFromDouble = eToWord32
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 18446744073709551615]@
instance Elevator Word64 where
  eToWord8 = dropDown
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown
  {-# INLINE eToWord16 #-}
  eToWord32 = dropDown
  {-# INLINE eToWord32 #-}
  eToWord64 = id
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1
  {-# INLINE eToDouble #-}
  eFromDouble = eToWord64
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 18446744073709551615]@ on 64bit
instance Elevator Word where
  eToWord8 = dropDown
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown
  {-# INLINE eToWord16 #-}
  eToWord32 = dropDown
  {-# INLINE eToWord32 #-}
  eToWord64 = fromIntegral
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 127]@
instance Elevator Int8 where
  eToWord8 = fromIntegral . max 0
  {-# INLINE eToWord8 #-}
  eToWord16 = raiseUp . max 0
  {-# INLINE eToWord16 #-}
  eToWord32 = raiseUp . max 0
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp . max 0
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1 . max 0
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1 . max 0
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 32767]@
instance Elevator Int16 where
  eToWord8 = dropDown . max 0
  {-# INLINE eToWord8 #-}
  eToWord16 = fromIntegral . max 0
  {-# INLINE eToWord16 #-}
  eToWord32 = raiseUp . max 0
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp . max 0
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1 . max 0
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1 . max 0
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 2147483647]@
instance Elevator Int32 where
  eToWord8 = dropDown . max 0
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown . max 0
  {-# INLINE eToWord16 #-}
  eToWord32 = fromIntegral . max 0
  {-# INLINE eToWord32 #-}
  eToWord64 = raiseUp . max 0
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1 . max 0
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1 . max 0
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 9223372036854775807]@
instance Elevator Int64 where
  eToWord8 = dropDown . max 0
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown . max 0
  {-# INLINE eToWord16 #-}
  eToWord32 = dropDown . max 0
  {-# INLINE eToWord32 #-}
  eToWord64 = fromIntegral . max 0
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1 . max 0
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1 . max 0
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0, 9223372036854775807]@ on 64bit
instance Elevator Int where
  eToWord8 = dropDown . max 0
  {-# INLINE eToWord8 #-}
  eToWord16 = dropDown . max 0
  {-# INLINE eToWord16 #-}
  eToWord32 = dropDown . max 0
  {-# INLINE eToWord32 #-}
  eToWord64 = fromIntegral . max 0
  {-# INLINE eToWord64 #-}
  eToFloat = squashTo1 . max 0
  {-# INLINE eToFloat #-}
  eToDouble = squashTo1 . max 0
  {-# INLINE eToDouble #-}
  eFromDouble = stretch . clamp01
  {-# INLINE eFromDouble #-}

-- | Values between @[0.0, 1.0]@
instance Elevator Float where
  eToWord8 = stretch . clamp01
  {-# INLINE eToWord8 #-}
  eToWord16 = stretch . clamp01
  {-# INLINE eToWord16 #-}
  eToWord32 = stretch . clamp01
  {-# INLINE eToWord32 #-}
  eToWord64 = stretch . clamp01
  {-# INLINE eToWord64 #-}
  eToFloat = id
  {-# INLINE eToFloat #-}
  eToDouble = float2Double
  {-# INLINE eToDouble #-}
  eFromDouble = eToFloat
  {-# INLINE eFromDouble #-}

-- | Values between @[0.0, 1.0]@
instance Elevator Double where
  eToWord8 = stretch . clamp01
  {-# INLINE eToWord8 #-}
  eToWord16 = stretch . clamp01
  {-# INLINE eToWord16 #-}
  eToWord32 = stretch . clamp01
  {-# INLINE eToWord32 #-}
  eToWord64 = stretch . clamp01
  {-# INLINE eToWord64 #-}
  eToFloat = double2Float
  {-# INLINE eToFloat #-}
  eToDouble = id
  {-# INLINE eToDouble #-}
  eFromDouble = id
  {-# INLINE eFromDouble #-}

-- | Discards imaginary part and changes precision of real part.
instance (Num e, Elevator e, RealFloat e) => Elevator (C.Complex e) where
  eToWord8 = eToWord8 . C.realPart
  {-# INLINE eToWord8 #-}
  eToWord16 = eToWord16 . C.realPart
  {-# INLINE eToWord16 #-}
  eToWord32 = eToWord32 . C.realPart
  {-# INLINE eToWord32 #-}
  eToWord64 = eToWord64 . C.realPart
  {-# INLINE eToWord64 #-}
  eToFloat = eToFloat . C.realPart
  {-# INLINE eToFloat #-}
  eToDouble = eToDouble . C.realPart
  {-# INLINE eToDouble #-}
  eFromDouble = (C.:+ 0) . eFromDouble
  {-# INLINE eFromDouble #-}

max_stars_repo_name: masterdezign/massiv

max_issues_repo_name: masterdezign/massiv

max_forks_repo_name: masterdezign/massiv


text:
import qualified System.Environment as SE
import qualified Data.ByteString.Lazy.Char8 as LC
import qualified Data.ByteString.Char8 as BC
import qualified Data.List as DL
import qualified Data.Map as DM
import qualified Data.Vector as DV
import qualified Statistics.Sample as SS
import qualified Statistics.Quantile as SQ
import Text.Printf

main :: IO()
main = do
    input <- LC.getContents
    output (nums $ LC.lines input)

-- clean input into a list of numbers
nums :: [LC.ByteString] -> [Double]
nums llns = floats $ map BC.unpack $ slns llns
    where
        slns = map (BC.concat . LC.toChunks)
        floats [] = []
        floats [x] = [read x :: Double]
        floats (x:xs) = [read x :: Double] ++ floats xs

-- | calculate the frequency of items in a set
-- [ [(number, frequency), ....
freqs xs = DM.toList $ DM.fromListWith (+) [(c, 1) | c <- xs]

-- | calculate the mode of a list of numbers
mode xs = (fst . DL.head) $ DL.sortBy sortGT $ freqs xs
    where
        sortGT (a1, b1) (a2, b2)
            | b1 < b2 = GT
            | b1 > b2 = LT
            | b1 == b2 = compare b1 b2

-- | calculate a quantile point
qtile q t xs = (SQ.continuousBy SQ.medianUnbiased q t) $ DV.fromList xs

-- | calculate the midhinge
midhinge xs = (qtile 1 4 xs + qtile 3 4 xs) / 2

-- | calculate the trimean
trimean xs = (qtile 2 4 xs + midhinge xs) / 2

h :: (Floating a, Ord t) => [t] -> a
-- | calculate Shannon Entropy
h xs = negate . sum . map (\(i,j) -> (p j xs) * (logBase 2 $ p j xs)) $ freqs xs
    where
        p n lst = (/) n tot_len
	tot_len = (fromIntegral $ length xs)

-- | string formatter: round to 4 decimal places
format f n = printf "%.4f" (f $ DV.fromList n)

-- | display output formatting
display xs = sequence_ [putStrLn (a++" : "++b) | (a,b) <- xs]

output :: [Double] -> IO ()
output xs = display table
    where table = [("Length", show $ length xs)
                  ,("Min", show $ minimum xs)
                  ,("Max", show $ maximum xs)
                  ,("Range", printf "%.4f" $ maximum xs - minimum xs)
                  ,("Q1", format (SQ.continuousBy SQ.medianUnbiased 1 4) xs)
                  ,("Q2", format (SQ.continuousBy SQ.medianUnbiased 2 4) xs)
                  ,("Q3", format (SQ.continuousBy SQ.medianUnbiased 3 4) xs)
                  ,("IQR", format (SQ.midspread SQ.medianUnbiased 4) xs)
                  ,("Trimean", printf "%.4f" $ trimean xs)
                  ,("Midhinge", printf "%.4f" $ midhinge xs)
                  ,("Mean", format SS.mean xs)
                  ,("Mode", show $ mode xs)
                  ,("Kurtosis", format SS.kurtosis xs)
                  ,("Skewness", format SS.skewness xs)
                  ,("Entropy", printf "%.4f" ((h xs) :: Float))]

max_stars_repo_name: drbunsen/describe

max_issues_repo_name: drbunsen/describe

max_forks_repo_name: drbunsen/describe


text:
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE OverloadedStrings #-}

module Main where

import Data.Monoid

import qualified Pipes.Prelude as P

import Statistics.Quantile.Util
import Statistics.Quantile.Exact
import Statistics.Quantile.Approximate.Sampling
import Statistics.Quantile.Bench
import Statistics.Quantile.Bench.Accuracy
import Statistics.Quantile.Types

import System.IO

data SelectorAccuracy = SelectorAccuracy String Deviation
  deriving (Eq, Show)

renderAccuracy :: SelectorAccuracy -> String
renderAccuracy (SelectorAccuracy s d) = s <> "," <> renderDeviation d

main :: IO ()
main = do
  hSetBuffering stdin LineBuffering
  hSetBuffering stdout LineBuffering
  (fp, fh) <- openTempFile tmpDir "accuracy.txt"
  n <- countOut fh (streamHandle stdin)
  hClose fh
  inh <- openFile fp ReadMode
  tv <- selectFromHandle median external inh
  hClose inh
  mapM (score fp tv) (estimators n sqrtN)
    >>= mapM_ (putStrLn . renderAccuracy)
  where estimators n f = [ ("exact", external)
                         , ("sampling", sampling (f n) n)
                         , ("jackknife-sampling", samplingJackknife (f n) n)
                         ]

        score fp tv (name, candidate) = do
          a <- benchAccuracy 10 fp median tv candidate
          pure (SelectorAccuracy name a)

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
{-# LANGUAGE OverloadedLists #-}
module Ch4 where

import Data.List
import Control.Monad
import Numeric.LinearAlgebra
import Graphics.Gnuplot.Simple

import Perceptron
import NeuralNetwork
import qualified MNIST
import SimpleNet (SimpleNet (..))
import qualified SimpleNet

main :: IO ()
main = do
  displayFunction1
  displayTangents
  displayFunction2
  displayGradients
  displayDescents

-- ch4_3_2
function1 :: Double -> Double
function1 x = 0.01 * x * x + 0.1 * x

displayFunction1 :: IO ()
displayFunction1 = plotPath attrs $ zip xs ys
  where attrs = [(Title "Function 1"), (YRange (0, 6))]
        xs = linearScale 100 (0, 20) :: [Double]
        ys = function1 <$> xs

tangent :: (Double -> Double) -> Double -> Double -> Double
tangent f x0 x = a * (x - x0) + f x0
  where a = numericalDiff f x0

displayTangents :: IO ()
displayTangents = plotPathsStyle attrs paths
  where attrs = [(Title "Tangents of Function 1"), (YRange (-1, 6))]
        xs = linearScale 100 (0, 20) :: [Double]
        paths = do
          (fn, name) <- [ (function1, "0.01 * x^2 + 0.1 * x")
                        , (tangent function1 5, "tangent (x = 5)")
                        , (tangent function1 10, "tangent (x = 10)")
                        ]
          let style = PlotStyle Lines $ CustomStyle [(LineTitle name)]
          return $ (style, zip xs (fn <$> xs))

-- ch4_3_3
function2 :: Vector R -> Double
function2 v = dot v v

displayFunction2 :: IO ()
displayFunction2 = plotFunc3d attrs [] xs ys function2'
  where
    attrs = [Title "Function 2"]
    xs = linearScale 100 (-3, 3) :: [Double]
    ys = linearScale 100 (-3, 3) :: [Double]
    function2' x y = function2 $ vector [x, y]

displayGradients :: IO ()
displayGradients = plotListStyle attrs style vs
  where
    attrs = [Title "Gradients of Function 2"]
    style = PlotStyle Vectors $ CustomStyle [(LineTitle "gradients")]
    xs = [(-2.0), (-1.75) .. 2.0] :: [Double]
    ys = [(-2.0), (-1.75) .. 2.0] :: [Double]
    vs = do
      x <- xs
      y <- ys
      let [dx, dy] = toList $ numericalGradient function2 $ vector [x, y]
      return ((x, y), (- dx / 20, - dy / 20))

displayDescents :: IO ()
displayDescents = plotPathStyle attrs style vs
  where
    attrs = [Title "Gradient descents of Function 2", XRange (-4, 4), YRange (-4, 4)]
    style = PlotStyle Points $ CustomStyle [(LineTitle "descents")]
    v0 = vector [-3.0, 4.0]
    vs = do
      v <- take 100 $ gradientDescents function2 0.1 v0
      let [x, y] = toList v
      return (x, y)

-- ch4_4_2

functionLoss :: Vector R -> Vector R -> Matrix R -> Double
functionLoss input label w = SimpleNet.loss net' input label
  where net' = SimpleNet w

net :: SimpleNet
net = SimpleNet $ (2 >< 3) [ 0.47355232, 0.99773930, 0.84668094
                           , 0.85557411, 0.03563661, 0.69422093]

input :: Vector R
input = [0.6, 0.9]              -- input

answer :: Vector R
answer = [0, 0, 1]              -- correct answer

functionNet :: Matrix R -> Double
functionNet = functionLoss input answer

dw :: SimpleNet -> Matrix R
dw (SimpleNet w) = numericalGradient functionNet w

netDescents :: SimpleNet -> [SimpleNet]
netDescents net = iterate descend net
  where
    lr = 0.1
    descend net@(SimpleNet w) = SimpleNet $ w - scale lr (dw net)

displayLearnings :: IO ()
displayLearnings = plotPathsStyle attrs paths
  where
    attrs = [(Title "Learnings of simple net"), (YRange (0, 1))]
    xs = [0 .. 50] :: [Double]
    n = length xs
    nets = netDescents net
    outputs = map (toList . softmax . (flip SimpleNet.predict input)) nets
    paths = do
      (i, path) <- zip [0..] $ transpose $ take n outputs
      let style = PlotStyle Lines $ CustomStyle [(LineTitle (show i))]
      return (style, zip xs path)

max_stars_repo_name: kayhide/deep-learning-from-scratch

max_issues_repo_name: kayhide/zero-to-deep-learning

max_forks_repo_name: kayhide/zero-to-deep-learning


text:
module GDAL.Complex (
    Complex(..)
  , module C
) where

import Data.Complex as C hiding (Complex(..))
import qualified Data.Complex as C

newtype Complex a =
  Complex { unComplex :: C.Complex a }
  deriving (Eq)

instance Show a => Show (Complex a) where
  show (Complex c) = show c

max_stars_repo_name: peddie/bindings-gdal

max_issues_repo_name: peddie/bindings-gdal

max_forks_repo_name: peddie/bindings-gdal


text:
{-# LANGUAGE CPP #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE DeriveTraversable #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE RankNTypes #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE DeriveDataTypeable #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}

#ifndef MIN_VERSION_hashable
#define MIN_VERSION_hashable(x,y,z) 1
#endif

-----------------------------------------------------------------------------
-- |
-- License     :  BSD-style (see the file LICENSE)
-- Maintainer  :  Edward Kmett <ekmett@gmail.com>
-- Stability   :  provisional
-- Portability :  portable
--
-- Operations on affine spaces.
-----------------------------------------------------------------------------
module Linear.Affine where

import Control.Applicative
import Control.DeepSeq
import Control.Monad (liftM)
import Control.Lens
import Data.Binary as Binary
import Data.Bytes.Serial
import Data.Coerce
import Data.Complex (Complex)
import Data.Data
import Data.Distributive
import Data.Foldable as Foldable
import Data.Functor.Bind
import Data.Functor.Classes
import Data.Functor.Product
import Data.Functor.Rep as Rep
import Data.HashMap.Lazy (HashMap)
import Data.Hashable
import Data.Hashable.Lifted
import Data.IntMap (IntMap)
import Data.Ix
import Data.Kind
import Data.Map (Map)
#if !(MIN_VERSION_base(4,11,0))
import Data.Semigroup (Semigroup)
#endif
import Data.Serialize as Cereal
import Data.Vector (Vector)
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed.Base as U
import Foreign.Storable
import GHC.Generics (Generic, Generic1)
import Linear.Epsilon
import Linear.Metric
import Linear.Plucker
import Linear.Quaternion
import Linear.V
import Linear.V0
import Linear.V1
import Linear.V2
import Linear.V3
import Linear.V4
import Linear.Vector
import System.Random (Random(..))

-- | An affine space is roughly a vector space in which we have
-- forgotten or at least pretend to have forgotten the origin.
--
-- > a .+^ (b .-. a)  =  b@
-- > (a .+^ u) .+^ v  =  a .+^ (u ^+^ v)@
-- > (a .-. b) ^+^ v  =  (a .+^ v) .-. q@
class Additive (Diff p) => Affine p where
  type Diff p :: Type -> Type

  infixl 6 .-.
  -- | Get the difference between two points as a vector offset.
  (.-.) :: Num a => p a -> p a -> Diff p a

  infixl 6 .+^
  -- | Add a vector offset to a point.
  (.+^) :: Num a => p a -> Diff p a -> p a

  infixl 6 .-^
  -- | Subtract a vector offset from a point.
  (.-^) :: Num a => p a -> Diff p a -> p a
  p .-^ v = p .+^ negated v
  {-# INLINE (.-^) #-}

instance (Affine f, Affine g) => Affine (Product f g) where
  type Diff (Product f g) = Product (Diff f) (Diff g)
  Pair a b .-. Pair c d = Pair (a .-. c) (b .-. d)
  Pair a b .+^ Pair c d = Pair (a .+^ c) (b .+^ d)
  Pair a b .-^ Pair c d = Pair (a .+^ c) (b .+^ d)

-- | Compute the quadrance of the difference (the square of the distance)
qdA :: (Affine p, Foldable (Diff p), Num a) => p a -> p a -> a
qdA a b = Foldable.sum (fmap (join (*)) (a .-. b))
{-# INLINE qdA #-}

-- | Distance between two points in an affine space
distanceA :: (Floating a, Foldable (Diff p), Affine p) => p a -> p a -> a
distanceA a b = sqrt (qdA a b)
{-# INLINE distanceA #-}

#define ADDITIVEC(CTX,T) instance CTX => Affine T where type Diff T = T ; \
  (.-.) = (^-^) ; {-# INLINE (.-.) #-} ; (.+^) = (^+^) ; {-# INLINE (.+^) #-} ; \
  (.-^) = (^-^) ; {-# INLINE (.-^) #-}
#define ADDITIVE(T) ADDITIVEC((), T)

ADDITIVE([])
ADDITIVE(Complex)
ADDITIVE(ZipList)
ADDITIVE(Maybe)
ADDITIVE(IntMap)
ADDITIVE(Identity)
ADDITIVE(Vector)
ADDITIVE(V0)
ADDITIVE(V1)
ADDITIVE(V2)
ADDITIVE(V3)
ADDITIVE(V4)
ADDITIVE(Plucker)
ADDITIVE(Quaternion)
ADDITIVE(((->) b))
ADDITIVEC(Ord k, (Map k))
ADDITIVEC((Eq k, Hashable k), (HashMap k))
ADDITIVEC(Dim n, (V n))

-- | A handy wrapper to help distinguish points from vectors at the
-- type level
newtype Point f a = P (f a)
  deriving ( Eq, Ord, Show, Read, Monad, Functor, Applicative, Foldable
           , Eq1, Ord1, Show1, Read1
           , Traversable, Apply, Additive, Metric
           , Fractional , Num, Ix, Storable, Epsilon
           , Semigroup, Monoid
           , Random, Hashable
           , Generic, Generic1, Data
           )

instance Finite f => Finite (Point f) where
  type Size (Point f) = Size f
  toV (P v) = toV v
  fromV v = P (fromV v)

instance NFData (f a) => NFData (Point f a) where
  rnf (P x) = rnf x

instance Serial1 f => Serial1 (Point f) where
  serializeWith f (P p) = serializeWith f p
  deserializeWith m = P `liftM` deserializeWith m

instance Serial (f a) => Serial (Point f a) where
  serialize (P p) = serialize p
  deserialize = P `liftM` deserialize

instance Binary (f a) => Binary (Point f a) where
  put (P p) = Binary.put p
  get = P `liftM` Binary.get

instance Serialize (f a) => Serialize (Point f a) where
  put (P p) = Cereal.put p
  get = P `liftM` Cereal.get

instance Hashable1 f => Hashable1 (Point f) where
  liftHashWithSalt h s (P f) = liftHashWithSalt h s f
  {-# INLINE liftHashWithSalt #-}

lensP :: Lens' (Point g a) (g a)
lensP afb (P a) = P <$> afb a
{-# INLINE lensP #-}

_Point :: Iso' (Point f a) (f a)
_Point = iso (\(P a) -> a) P
{-# INLINE _Point #-}

instance (t ~ Point g b) => Rewrapped (Point f a) t
instance Wrapped (Point f a) where
  type Unwrapped (Point f a) = f a
  _Wrapped' = _Point
  {-# INLINE _Wrapped' #-}

-- These are stolen from Data.Profunctor.Unsafe
(.#) :: Coercible b a => (b -> c) -> (a -> b) -> a -> c
f .# _ = coerce f
{-# INLINE (.#) #-}

(#.) :: Coercible c b => (b -> c) -> (a -> b) -> a -> c
(#.) _ = coerce (\x -> x :: b) :: forall a b. Coercible b a => a -> b
{-# INLINE (#.) #-}

unP :: Point f a -> f a
unP (P x) = x
{-# INLINE unP #-}

-- We can't use GND to derive 'Bind' because 'join' causes
-- role troubles. However, GHC 7.8 and above let us use
-- explicit coercions for (>>-).
instance Bind f => Bind (Point f) where
  (>>-) = ((P .) . (. (unP .))) #. (>>-) .# unP
  join (P m) = P $ m >>- \(P m') -> m'

instance Distributive f => Distributive (Point f) where
  distribute = P . collect (\(P p) -> p)
  collect = (P .) #. collect .# (unP .)

instance Representable f => Representable (Point f) where
  type Rep (Point f) = Rep f
  tabulate = P #. tabulate
  {-# INLINE tabulate #-}
  index = Rep.index .# unP
  {-# INLINE index #-}

type instance Index (Point f a) = Index (f a)
type instance IxValue (Point f a) = IxValue (f a)

instance Ixed (f a) => Ixed (Point f a) where
  ix l = lensP . ix l
  {-# INLINE ix #-}

instance Traversable f => Each (Point f a) (Point f b) a b where
  each = traverse
  {-# INLINE each #-}

instance R1 f => R1 (Point f) where
  _x = lensP . _x
  {-# INLINE _x #-}

instance R2 f => R2 (Point f) where
  _y = lensP . _y
  {-# INLINE _y #-}
  _xy = lensP . _xy
  {-# INLINE _xy #-}

instance R3 f => R3 (Point f) where
  _z = lensP . _z
  {-# INLINE _z #-}
  _xyz = lensP . _xyz
  {-# INLINE _xyz #-}

instance R4 f => R4 (Point f) where
  _w = lensP . _w
  {-# INLINE _w #-}
  _xyzw = lensP . _xyzw
  {-# INLINE _xyzw #-}

instance Additive f => Affine (Point f) where
  type Diff (Point f) = f
  (.-.) = (. unP) #. (^-^) .# unP
  {-# INLINE (.-.) #-}
  (.+^) = (P .) #. (^+^) .# unP
  {-# INLINE (.+^) #-}
  (.-^) = (P .) #. (^-^) .# unP
  {-# INLINE (.-^) #-}

-- | Vector spaces have origins.
origin :: (Additive f, Num a) => Point f a
origin = P zero

-- | An isomorphism between points and vectors, given a reference
--   point.
relative :: (Additive f, Num a) => Point f a -> Iso' (Point f a) (f a)
relative p0 = iso (.-. p0) (p0 .+^)
{-# INLINE relative #-}

newtype instance U.Vector    (Point f a) =  V_P (U.Vector    (f a))
newtype instance U.MVector s (Point f a) = MV_P (U.MVector s (f a))
instance U.Unbox (f a) => U.Unbox (Point f a)

instance U.Unbox (f a) => M.MVector U.MVector (Point f a) where
  {-# INLINE basicLength #-}
  {-# INLINE basicUnsafeSlice #-}
  {-# INLINE basicOverlaps #-}
  {-# INLINE basicUnsafeNew #-}
  {-# INLINE basicUnsafeRead #-}
  {-# INLINE basicUnsafeWrite #-}
  basicLength (MV_P v) = M.basicLength v
  basicUnsafeSlice m n (MV_P v) = MV_P (M.basicUnsafeSlice m n v)
  basicOverlaps (MV_P v) (MV_P u) = M.basicOverlaps v u
  basicUnsafeNew n = MV_P `liftM` M.basicUnsafeNew n
  basicUnsafeRead (MV_P v) i = P `liftM` M.basicUnsafeRead v i
  basicUnsafeWrite (MV_P v) i (P x) = M.basicUnsafeWrite v i x
  basicInitialize (MV_P v) = M.basicInitialize v
  {-# INLINE basicInitialize #-}

instance U.Unbox (f a) => G.Vector U.Vector (Point f a) where
  {-# INLINE basicUnsafeFreeze #-}
  {-# INLINE basicUnsafeThaw   #-}
  {-# INLINE basicLength       #-}
  {-# INLINE basicUnsafeSlice  #-}
  {-# INLINE basicUnsafeIndexM #-}
  basicUnsafeFreeze (MV_P v) = V_P `liftM` G.basicUnsafeFreeze v
  basicUnsafeThaw   ( V_P v) = MV_P `liftM` G.basicUnsafeThaw   v
  basicLength       ( V_P v) = G.basicLength v
  basicUnsafeSlice m n (V_P v) = V_P (G.basicUnsafeSlice m n v)
  basicUnsafeIndexM (V_P v) i = P `liftM` G.basicUnsafeIndexM v i

max_stars_repo_name: ekmett/linear

max_issues_repo_name: ekmett/linear

max_forks_repo_name: ekmett/linear