text:
-----------------------------------------------------------------------------
--
-- Module      :  Drool.Utils.RenderHelpers
-- Copyright   :  Tobias Fuchs
-- License     :  AllRightsReserved
--
-- Maintainer  :  twh.fuchs@gmail.com
-- Stability   :  experimental
-- Portability :  POSIX
--
-- |
--
-----------------------------------------------------------------------------

{-# OPTIONS -O2 -Wall #-}

module Drool.Utils.RenderHelpers (
    RenderSettings(..),
    applyGlobalRotation,
    applyPerspective,
    nextPerspective,
    useLight,
    applyBandRangeAmp,
    bandRangeAmpSamples,
    scaleSamples,

    vertexWithNormal,
    vertexWithNormalAndColor,
    normalsFromVertices,
    drawNormal,
    color3AddAlpha,
    color4MulAlpha,
    color3MulValue,
    color4MulValue,

    vertexToVector,
    vx4ToVx3,

    v3x,
    v3y,
    v3z,
    vx3x,
    vx3y,
    vx3z,
    vx4x,
    vx4y,
    vx4z,
    vx4w,

    getViewpointFromModelView,
    featuresToIntensity
) where

import Debug.Trace

-- Imports
-- {{{
import Data.IORef ( IORef )

import Data.Packed.Matrix as HM ( (><), Matrix )
import Numeric.LinearAlgebra.Algorithms as LA ( inv )
import Numeric.Container as NC ( vXm, fromList, toList )

import qualified Drool.Types as DT ( SignalList(..), RenderPerspective(..), RotationVector(..) )
import Drool.Utils.SigGen as SigGen ( SValue, TValue, SignalGenerator(..) )
import Drool.Utils.FeatureExtraction as FE (
    SignalFeatures(..), SignalFeaturesList(..),
    emptyFeatures,
    FeatureTarget(..), featureTargetFromIndex )
import Drool.Utils.Conversions as Conv ( interleaveArrays, aZip )
import Drool.ApplicationContext as AC ( ContextSettings(..), LightConfig(..) )
import Graphics.Rendering.OpenGL (
    Vector3 (..),
    Vertex3 (..),
    Vertex4 (..),
    Normal3 (..),
    Color3(..),
    Color4(..),
    Light(..),
    light,
    VertexComponent,
    renderPrimitive,
    PrimitiveMode(..),
    vertex, normal,
    color,
    materialEmission,
    materialAmbient,
    materialDiffuse,
    materialSpecular,
    materialShininess,
    translate,
    MatrixOrder(..),
    getMatrixComponents,
    GLmatrix,
    ($=),
    Face(..),
    Capability(..),
    polygonOffsetFill,
    polygonOffsetLine,
    ambient,
    diffuse,
    specular,
    rotate,
    translate,
    GLfloat )
import qualified Graphics.Rendering.FTGL as FTGL
import qualified Control.Concurrent.MVar as MV ( MVar )
import qualified Control.Concurrent.Chan as CC
-- }}}

-- Settings rendering independent from visual used:
data RenderSettings = RenderSettings { signalGenerator :: SignalGenerator,
                                       samplingSem :: MV.MVar Int,
                                       renderingSem :: MV.MVar Int,
                                       numNewSignalsChan :: CC.Chan Int,
                                       -- IORef to signal buffer:
                                       signalBuf :: IORef DT.SignalList,
                                       -- Position of light 0
                                       lightPos0 :: Vertex4 GLfloat,
                                       -- Position of light 1
                                       lightPos1 :: Vertex4 GLfloat,
                                       -- Containing one SignalFeatures component for every signal:
                                       featuresBuf :: IORef (FE.SignalFeaturesList),
                                       -- Current number of signals in signal buffer:
                                       numSignals :: Int,
                                       -- Number of samples in most recent signal:
                                       numSamples :: Int,
                                       -- Number of new signals since last render pass:
                                       numNewSignals :: Int,
                                       -- Whether to reverse the signal buffer:
                                       reverseBuffer :: Bool,
                                       -- Tick of rendering pass:
                                       tick :: Int }

applyPerspective :: DT.RenderPerspective -> IO ()
applyPerspective p = do
  case p of
    DT.Isometric -> do
      rotate (45::GLfloat) $ Vector3 1.0 0.0 0.0
      rotate (45::GLfloat) $ Vector3 0.0 1.0 0.0
    DT.Top -> do
      rotate (90::GLfloat) $ Vector3 1.0 0.0 0.0
    DT.Front -> do
      rotate (20.0::GLfloat) $ Vector3 1.0 0.0 0.0
    DT.Side -> do
      rotate (20.0::GLfloat) $ Vector3 1.0 0.0 0.0
      rotate (-90::GLfloat) $ Vector3 0.0 1.0 0.0

applyGlobalRotation :: DT.RotationVector -> DT.RotationVector -> IO ()
applyGlobalRotation fixedRotation incRotation = do
  rotate (DT.rotX fixedRotation) $ Vector3 1.0 0.0 0.0
  rotate (DT.rotX incRotation)   $ Vector3 1.0 0.0 0.0
  rotate (DT.rotY fixedRotation) $ Vector3 0.0 1.0 0.0
  rotate (DT.rotY incRotation)   $ Vector3 0.0 1.0 0.0
  rotate (DT.rotZ fixedRotation) $ Vector3 0.0 0.0 1.0
  rotate (DT.rotZ incRotation)   $ Vector3 0.0 0.0 1.0

nextPerspective :: DT.RenderPerspective -> DT.RenderPerspective
nextPerspective cur = case cur of
  DT.Isometric -> DT.Top
  DT.Top       -> DT.Front
  DT.Front     -> DT.Side
  DT.Side      -> DT.Isometric

useLight :: AC.LightConfig -> IO ()
useLight lightConfig = do
  let lightState = AC.lightState lightConfig
  let lightIdx   = fromIntegral $ AC.lightIndex lightConfig
  light (Light lightIdx) $= lightState
  let lighting state =
        if state == Enabled then do
          let intensity = AC.lightIntensity lightConfig
          ambient  (Light lightIdx) $= color4MulValue (AC.lightAmbient  lightConfig) intensity
          diffuse  (Light lightIdx) $= color4MulValue (AC.lightDiffuse  lightConfig) intensity
          specular (Light lightIdx) $= color4MulValue (AC.lightSpecular lightConfig) intensity
        else
          return ()
  lighting lightState

-- Expects a sample, t, number of samples in total, list of band range amplifiers, and returns amplified sample for t.
applyBandRangeAmp :: SValue -> TValue -> Int -> [Float] -> SValue
-- {{{
applyBandRangeAmp s t nSamples amps = s * ampValue
  where numRanges      = length amps
        rangeWidth     = nSamples `div` numRanges                            -- ...[-------]...
        rangePos       = t `mod` rangeWidth                                  -- ...|....x..|...
        activeRangeIdx = max 0 $ ((fromIntegral t)-rangePos) `div` rangeWidth :: Int -- [.|x|.|.|.]
        -- Amp value of active range
        amp_1          = if length amps > activeRangeIdx then realToFrac $ amps !! activeRangeIdx else 1.0
        -- Amp value of next range, use range n as next range of range n
        nextRangeIdx   = max 0 $ min (activeRangeIdx+1) (numRanges-1)
        amp_2          = if length amps > nextRangeIdx then realToFrac $ amps !! nextRangeIdx else 1.0
        ampValue       = amp_1 + ((amp_2-amp_1) * (fromIntegral (rangePos) / fromIntegral (rangeWidth-1)))
-- }}}

bandRangeAmpSamples :: [SValue] -> [Float] -> [SValue]
bandRangeAmpSamples samples amps = bandRangeAmpSamplesRec samples amps (length samples) 0

bandRangeAmpSamplesRec :: [SValue] -> [Float] -> Int -> TValue -> [SValue]
bandRangeAmpSamplesRec (x:xs) amps nSamples t = ampSample : (bandRangeAmpSamplesRec xs amps nSamples (t+1))
  where ampSample = (applyBandRangeAmp x t nSamples amps)
bandRangeAmpSamplesRec [] _ _ _ = []

scaleSamples :: (Fractional a) => [a] -> a -> [a]
scaleSamples samples a = map ( \s -> s * a ) samples

-- {{{

-- Resolve x component of a 3-dimensional Vector
v3x :: Vector3 a -> a
v3x (Vector3 x _ _) = x
-- Resolve y component of a 3-dimensional Vector
v3y :: Vector3 a -> a
v3y (Vector3 _ y _) = y
-- Resolve z component of a 3-dimensional Vector
v3z :: Vector3 a -> a
v3z (Vector3 _ _ z) = z

-- Resolve x component of a 3-dimensional Vertex
vx3x :: Vertex3 a -> a
vx3x (Vertex3 x _ _) = x
-- Resolve y component of a 3-dimensional Vertex
vx3y :: Vertex3 a -> a
vx3y (Vertex3 _ y _) = y
-- Resolve z component of a 3-dimensional Vertex
vx3z :: Vertex3 a -> a
vx3z (Vertex3 _ _ z) = z

-- Resolve x component of a 4-dimensional Vertex
vx4x :: Vertex4 a -> a
vx4x (Vertex4 x _ _ _) = x
-- Resolve y component of a 4-dimensional Vertex
vx4y :: Vertex4 a -> a
vx4y (Vertex4 _ y _ _) = y
-- Resolve z component of a 4-dimensional Vertex
vx4z :: Vertex4 a -> a
vx4z (Vertex4 _ _ z _) = z
-- Resolve w component of a 4-dimensional Vertex
vx4w :: Vertex4 a -> a
vx4w (Vertex4 _ _ _ w) = w

-- Resolve x component of a 3-dimensional Normal
n3x :: Normal3 a -> a
n3x (Normal3 x _ _) = x
-- Resolve y component of a 3-dimensional Normal
n3y :: Normal3 a -> a
n3y (Normal3 _ y _) = y
-- Resolve z component of a 3-dimensional Normal
n3z :: Normal3 a -> a
n3z (Normal3 _ _ z) = z

-- TODO: All these vector / matrix operations should be
--       outsourced to BLAS.

vx4ToVx3 :: Vertex4 a -> Vertex3 a
vx4ToVx3 (Vertex4 x y z _) = Vertex3 x y z

vertexToVector :: (Num a) => Vertex3 a -> Vector3 a
vertexToVector v = Vector3 (vx3x v) (vx3y v) (vx3z v)

-- Cross product of two vectors
v3cross :: (Num a) => Vector3 a -> Vector3 a -> Vector3 a
v3cross (Vector3 a1 a2 a3) (Vector3 b1 b2 b3) = Vector3 c1 c2 c3
  where c1 = a2 * b3 - a3 * b2
        c2 = a3 * b1 - a1 * b3
        c3 = a1 * b2 - a2 * b1

v3add :: (Num a) => Vector3 a -> Vector3 a -> Vector3 a
v3add a b = Vector3 x y z
  where x = v3x a + v3x b
        y = v3y a + v3y b
        z = v3z a + v3z b
v3sub :: (Num a) => Vector3 a -> Vector3 a -> Vector3 a
v3sub a b = Vector3 x y z
  where x = v3x a - v3x b
        y = v3y a - v3y b
        z = v3z a - v3z b
v3sum :: (Num a) => [Vector3 a] -> Vector3 a
v3sum vlist = foldl (\accum v -> v3add v accum) (Vector3 0 0 0) vlist

v3div :: (Fractional a, Num a) => Vector3 a -> Vector3 a -> Vector3 a
v3div a b = Vector3 x y z
  where x = v3x a / v3x b
        y = v3y a / v3y b
        z = v3z a / v3z b

vx3add :: (Num a) => Vertex3 a -> Vertex3 a -> Vertex3 a
vx3add a b = Vertex3 x y z
  where x = vx3x a + vx3x b
        y = vx3y a + vx3y b
        z = vx3z a + vx3z b

n3Invert :: Normal3 GLfloat -> Normal3 GLfloat
n3Invert n = Normal3 x y z
  where x = -(n3x n)
        y = -(n3y n)
        z = -(n3z n)

color3AddAlpha :: Color3 GLfloat -> GLfloat -> Color4 GLfloat
color3AddAlpha (Color3 r g b) a = Color4 r g b a

color4MulAlpha :: Color4 GLfloat -> GLfloat -> Color4 GLfloat
color4MulAlpha (Color4 r g b a) x = Color4 r g b (a*x)

color4MulValue :: Color4 GLfloat -> GLfloat -> Color4 GLfloat
color4MulValue (Color4 r g b a) x = Color4 (r*x) (g*x) (b*x) a

color3MulValue :: Color3 GLfloat -> GLfloat -> Color3 GLfloat
color3MulValue (Color3 r g b) x = Color3 (r*x) (g*x) (b*x)

-- }}}

-- Useful for mapping over a list containing tuples of (vertex, vertexNormal):
vertexWithNormal :: (Vertex3 GLfloat, Normal3 GLfloat) -> IO ()
vertexWithNormal (v,n) = do normal n
                            vertex v
                            return ()

vertexWithNormalAndColor :: (Vertex3 GLfloat, Normal3 GLfloat) -> Color4 GLfloat -> IO ()
vertexWithNormalAndColor (v,n) c = do color c
                                      normal n
                                      vertex v
                                      return ()

drawNormal :: (Fractional a, VertexComponent a) => ( Vertex3 a, Normal3 a ) -> IO ()
drawNormal (v,n) = do let from = v
                      let nx   = (n3x n) * 100.0
                      let ny   = (n3y n) * 100.0
                      let nz   = (n3z n) * 100.0
                      let to   = Vertex3 (nx + vx3x v) (ny + vx3y v) (nz + vx3z v)
                      color $ (Color4 1.0 1.0 1.0 1.0 :: Color4 GLfloat)
                      renderPrimitive Lines ( do vertex from
                                                 vertex to )

normalsFromVertices :: [[ Vertex3 GLfloat ]] -> [ Normal3 GLfloat ]
-- {{{
normalsFromVertices sigs = normalsFromVertices' sigs 0

normalsFromVertices' :: [[ Vertex3 GLfloat ]] -> Int -> [ Normal3 GLfloat ]
normalsFromVertices' sigs xIdx = case sigs of
  (sigPrev:sig:sigNext:[]) -> if nSamples > 0 && xIdx < nSamples then resultNormal : normalsFromVertices' sigs (xIdx+1) else []
    where nSamples = length sig
          boundx   = \a -> max 0 (min a (nSamples-1)) :: Int
          vertexC  = sig !! xIdx
          vertexR  = sig !! (boundx (xIdx+1))
          vertexL  = sig !! (boundx (xIdx-1))
          vertexT  = if length sigPrev > xIdx then sigPrev !! xIdx else vertexC
          vertexB  = if length sigNext > xIdx then sigNext !! xIdx else vertexC
          point    = vertexToVector vertexC
          pointR   = vertexToVector vertexR
          pointL   = vertexToVector vertexL
          pointT   = vertexToVector vertexT
          pointB   = vertexToVector vertexB
          vR = v3sub pointR point   --   n4  T  n1
          vT = v3sub pointT point   --       |
          vL = v3sub pointL point   --    L--C--R
          vB = v3sub pointB point   --       |
          n1 = v3cross vT vR        --   n3  B  n2
          n2 = v3cross vR vB
          n3 = v3cross vB vL
          n4 = v3cross vL vT
          -- no normalization here as GL.normalize is enabled
          n  = v3div (v3sum [ n1, n2, n3, n4 ]) (Vector3 4.0 4.0 (4.0 :: GLfloat))
          resultNormal = n3Invert $ Normal3 (v3x n) (v3y n) (v3z n)
  _ -> []
-- }}}

getViewpointFromModelView :: GLmatrix GLfloat -> IO ( Vector3 GLfloat )
-- {{{
getViewpointFromModelView mvMatrix = do
  mvMatrixRows <- getMatrixComponents ColumnMajor mvMatrix
  -- OpenGL float matrix to GSL double matrix:
  let mvMatrixGSL = (4 >< 4) $ map (\e -> realToFrac e :: Double) mvMatrixRows
  -- Invert GSL matrix:
  let mvMatrixGSLInv = LA.inv ( mvMatrixGSL :: HM.Matrix Double )
  -- 4th column vector is position of view point. First 3 columns are
  -- transformation elements.
  let viewPointProjection = NC.fromList [ 0.0, 0.0, 0.0, 1.0 ]
  let viewPointModelView = NC.toList $ NC.vXm viewPointProjection mvMatrixGSLInv
  let viewPoint = Vector3 (realToFrac $ viewPointModelView !! 0)
                          (realToFrac $ viewPointModelView !! 1)
                          (realToFrac $ viewPointModelView !! 2) :: Vector3 GLfloat
  return viewPoint
-- }}}

featuresToIntensity :: FE.SignalFeatures -> FE.FeatureTarget -> AC.ContextSettings -> (GLfloat,GLfloat)
featuresToIntensity features target cSettings = (lCoeff, bCoeff)
  where loudness     = realToFrac $ FE.totalEnergy features
        basslevel    = realToFrac $ FE.bassEnergy features
        lTarget      = FE.featureTargetFromIndex $ AC.featureSignalEnergyTargetIdx cSettings
        bTarget      = FE.featureTargetFromIndex $ AC.featureBassEnergyTargetIdx cSettings
        lCoeff       = if lTarget == target || lTarget == FE.GlobalAndLocalTarget then (
                          realToFrac $ (AC.featureSignalEnergySurfaceCoeff cSettings) * loudness )
                       else 0.0
        bCoeff       = if bTarget == target || bTarget == FE.GlobalAndLocalTarget then (
                          realToFrac $ (AC.featureBassEnergySurfaceCoeff cSettings) * basslevel )
                       else 0.0

max_stars_repo_name: fuchsto/drool

max_issues_repo_name: fuchsto/drool

max_forks_repo_name: fuchsto/drool


text:
{-|
  Safe immutable intermediate (functional) Matlab data structures.

  This provides an alternative representation of Matlab data structures which can be marshalled to and from 'MXArray'.
  An instance of 'Show' is provided which creates Matlab expressions.
-}
module Foreign.Matlab.Array.IMX (
    IMXData(..),
    IMXArray,
    IMXArrayElem (imxConstr, imxArray),

    -- * Interface with "Foreign.Matlab.Array"
    imxData, iMXData,
    imxFun, iMXFun,
    withIMXData, takeIMXData,
    withIMXDataList,

    -- * Construction and access
    imxSize,
    listIMX, scalarIMX,
    imxList, imxScalar,
    listIMXStruct,
    imxStructList
  ) where

import           Control.Monad
import           Data.Array.IArray
import           Data.Complex
import           Data.List
import qualified Data.Map.Strict as DM
import           Foreign.Matlab.Util
import           Foreign.Matlab.Internal
import           Foreign.Matlab.Types
import           Foreign.Matlab.Array

type IMXArray a = Array MIndex a
-- |The basic immutable (functional) representation of Matlab data structures, representing a generic 'MXArray'
data IMXData =
    IMXNull
  | IMXCell     (IMXArray IMXData)
  | IMXStruct   [String] (Array (MIndex,Int) IMXData) -- ^ field name list and array mapping (index,field index) to values
  | IMXLogical  (IMXArray MLogical)
  | IMXChar     (IMXArray MChar)
  | IMXDouble   (IMXArray MDouble)
  | IMXSingle   (IMXArray MSingle)
  | IMXInt8     (IMXArray MInt8)
  | IMXUint8    (IMXArray MUint8)
  | IMXInt16    (IMXArray MInt16)
  | IMXUint16   (IMXArray MUint16)
  | IMXInt32    (IMXArray MInt32)
  | IMXUint32   (IMXArray MUint32)
  | IMXInt64    (IMXArray MInt64)
  | IMXUint64   (IMXArray MUint64)
  | IMXComplexDouble    (IMXArray (MComplex MDouble))
  | IMXComplexSingle    (IMXArray (MComplex MSingle))
  | IMXObject   String IMXData -- ^ object class name and object data, currently always IMXStruct
  deriving (Eq)

type IMXFun = [IMXData] -> Int -> IO [IMXData]

fixMSize :: MSize -> [a] -> MSize
fixMSize s l = maybe s (\i -> replaceIndex s i (length l `div` negate (product s))) $ elemIndex (-1) s

-- |Create an sized array for 'IMXData' from a sequential list of elements.  'MSize' may contain at most one -1 value, which will be inferred from the length of the list.
listIMXArray :: MSize -> [e] -> IMXArray e
listIMXArray s l = listArray (mSizeRange (fixMSize s l)) l

-- |Create a scalar (array) for 'IMXData'
scalarIMXArray :: e -> IMXArray e
scalarIMXArray e = listIMXArray [] [e]

class IMXArrayElem a where
  -- |Generic 'IMXData' 'IMXArray' constructor
  imxConstr :: IMXArray a -> IMXData
  -- |Generic 'IMXData' 'IMXArray' accessor
  imxArray :: IMXData -> Maybe (IMXArray a)
instance IMXArrayElem IMXData  where { imxConstr = IMXCell    ; imxArray (IMXCell    a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MLogical where { imxConstr = IMXLogical ; imxArray (IMXLogical a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MChar    where { imxConstr = IMXChar    ; imxArray (IMXChar    a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MDouble  where { imxConstr = IMXDouble  ; imxArray (IMXDouble  a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MSingle  where { imxConstr = IMXSingle  ; imxArray (IMXSingle  a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MInt8    where { imxConstr = IMXInt8    ; imxArray (IMXInt8    a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MUint8   where { imxConstr = IMXUint8   ; imxArray (IMXUint8   a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MInt16   where { imxConstr = IMXInt16   ; imxArray (IMXInt16   a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MUint16  where { imxConstr = IMXUint16  ; imxArray (IMXUint16  a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MInt32   where { imxConstr = IMXInt32   ; imxArray (IMXInt32   a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MUint32  where { imxConstr = IMXUint32  ; imxArray (IMXUint32  a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MInt64   where { imxConstr = IMXInt64   ; imxArray (IMXInt64   a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem MUint64  where { imxConstr = IMXUint64  ; imxArray (IMXUint64  a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem (MComplex MDouble) where { imxConstr = IMXComplexDouble ; imxArray (IMXComplexDouble a) = Just a ; imxArray _ = Nothing }
instance IMXArrayElem (MComplex MSingle) where { imxConstr = IMXComplexSingle ; imxArray (IMXComplexSingle a) = Just a ; imxArray _ = Nothing }

-- |Generic 'IMXData' list constructor.  Specified 'MSize' may contain at most one -1 value, which will be inferred from the length of the list.
listIMX :: IMXArrayElem a => MSize -> [a] -> IMXData
listIMX s = imxConstr . listIMXArray s
-- |Generic 'IMXData' scalar constructor
scalarIMX :: IMXArrayElem a => a -> IMXData
scalarIMX = imxConstr . scalarIMXArray
--emptyIMX :: IMXArrayElem a => a -> IMXData
--emptyIMX a = listIMX [0] ([] `asTypeOf` [a])

-- |Generic 'IMXData' array size accessor
imxSize :: IMXData -> MSize
imxSize IMXNull = [0]
imxSize (IMXCell    a) = mRangeSize (bounds a)
imxSize (IMXLogical a) = mRangeSize (bounds a)
imxSize (IMXChar    a) = mRangeSize (bounds a)
imxSize (IMXDouble  a) = mRangeSize (bounds a)
imxSize (IMXSingle  a) = mRangeSize (bounds a)
imxSize (IMXInt8    a) = mRangeSize (bounds a)
imxSize (IMXUint8   a) = mRangeSize (bounds a)
imxSize (IMXInt16   a) = mRangeSize (bounds a)
imxSize (IMXUint16  a) = mRangeSize (bounds a)
imxSize (IMXInt32   a) = mRangeSize (bounds a)
imxSize (IMXUint32  a) = mRangeSize (bounds a)
imxSize (IMXInt64   a) = mRangeSize (bounds a)
imxSize (IMXUint64  a) = mRangeSize (bounds a)
imxSize (IMXComplexDouble a) = mRangeSize (bounds a)
imxSize (IMXComplexSingle a) = mRangeSize (bounds a)
imxSize (IMXStruct _ a) = mRangeSize (r0,r1) where ((r0,_),(r1,_)) = bounds a
imxSize (IMXObject _ d) = imxSize d

-- |Generic 'IMXData' list accessor
imxList :: IMXArrayElem a => IMXData -> Maybe [a]
imxList = fmap elems . imxArray
-- |Generic 'IMXData' scalar accessor
imxScalar :: IMXArrayElem a => IMXData -> Maybe a
imxScalar a = case imxList a of
  Just [x] -> Just x
  _ -> Nothing

-- |Create a sized struct array from a sequential list of consecutive field values (@[i0f0,i0f1,...,i0fM,i1f0,i1f1,....,iNfM]@).
listIMXStruct :: [String] -> MSize -> [IMXData] -> IMXData
listIMXStruct f s l = IMXStruct f $ listArray r l where
  n = length f
  r = ((r0,0),(r1,pred n))
  (r0,r1) = mSizeRange $ tail $ fixMSize (n:s) l

-- |Access a struct as list of fields and list of consecutive field values
imxStructList :: IMXData -> Maybe ([String], [IMXData])
imxStructList (IMXStruct f v) = Just (f, elems v)
imxStructList _ = Nothing

-- |Create an immutable representation from an 'MXArray'
imxData :: MXArray a -> IO IMXData
imxData a = do
  t <- mxArrayClass a
  c <- if t /= MXClassNull then mxArrayIsComplex a else return undefined
  imxc t c where

  imxc :: MXClass -> Bool -> IO IMXData
  imxc MXClassNull _            = return IMXNull
  imxc MXClassCell False        = IMXCell       =.< imxa (imxData . mCell)
  imxc MXClassStruct False      = do
    s <- mxArraySize a'
    fv <- mxArrayGetAll a'
    f <- if null fv then mStructFields a' else return (map fst (DM.toList $ _mStruct (head fv)))
    listIMXStruct f s =.< mapM imxData (concatMap (map snd . DM.toList . _mStruct) fv)
  imxc MXClassLogical False     = IMXLogical    =.< imxa return
  imxc MXClassChar False        = IMXChar       =.< imxa return
  imxc MXClassDouble False      = IMXDouble     =.< imxa return
  imxc MXClassSingle False      = IMXSingle     =.< imxa return
  imxc MXClassInt8 False        = IMXInt8       =.< imxa return
  imxc MXClassUint8 False       = IMXUint8      =.< imxa return
  imxc MXClassInt16 False       = IMXInt16      =.< imxa return
  imxc MXClassUint16 False      = IMXUint16     =.< imxa return
  imxc MXClassInt32 False       = IMXInt32      =.< imxa return
  imxc MXClassUint32 False      = IMXUint32     =.< imxa return
  imxc MXClassInt64 False       = IMXInt64      =.< imxa return
  imxc MXClassUint64 False      = IMXUint64     =.< imxa return
  imxc MXClassDouble True       = IMXComplexDouble      =.< imxa return
  imxc MXClassSingle True       = IMXComplexSingle      =.< imxa return
  imxc MXClassObject False      = do
    Just c <- mObjectGetClass a'
    IMXObject c =.< imxc MXClassStruct False
  imxc t c = fail ("imxData: unhandled mxArray type " ++ show t ++ if c then "(complex)" else "")

  imxa :: MXArrayComponent a => (a -> IO b) -> IO (IMXArray b)
  imxa f = do
    s <- mxArraySize a'
    listIMXArray s =.< mapM f =<< mxArrayGetAll a'

  a' = unsafeCastMXArray a

-- |Create a new 'MXArray' from a functional representation.
iMXData :: IMXData -> IO MAnyArray
iMXData = imxd where
  imxd :: IMXData -> IO MAnyArray
  imxd IMXNull = return $ anyMXArray mNullArray
  imxd (IMXCell a) = imxa a (MCell .=< iMXData)
  imxd (IMXStruct f a) = do
    let ((r0,_),(r1,_)) = bounds a
    m <- createStruct (mRangeSize (r0,r1)) f
    zipWithM_ (\i -> mStructSetFields m (mOffset i) <=< mapM iMXData) [0..] (segment (length f) (elems a))
    return $ anyMXArray m
  imxd (IMXLogical a)   = imxa a return
  imxd (IMXChar a)      = imxa a return
  imxd (IMXDouble a)    = imxa a return
  imxd (IMXSingle a)    = imxa a return
  imxd (IMXInt8 a)      = imxa a return
  imxd (IMXUint8 a)     = imxa a return
  imxd (IMXInt16 a)     = imxa a return
  imxd (IMXUint16 a)    = imxa a return
  imxd (IMXInt32 a)     = imxa a return
  imxd (IMXUint32 a)    = imxa a return
  imxd (IMXInt64 a)     = imxa a return
  imxd (IMXUint64 a)    = imxa a return
  imxd (IMXComplexDouble a) = imxa a return
  imxd (IMXComplexSingle a) = imxa a return
  imxd (IMXObject c a) = do
    m <- imxd a
    mObjectSetClass (unsafeCastMXArray m) c
    return m

  imxa :: MXArrayComponent b => IMXArray a -> (a -> IO b) -> IO MAnyArray
  imxa a f = do
    m <- createMXArray (mRangeSize (bounds a))
    mxArraySetAll m =<< mapM f (elems a)
    return $ anyMXArray m

withIMXData :: IMXData -> (MAnyArray -> IO a) -> IO a
withIMXData d f = do
  a <- iMXData d
  r <- f a
  freeMXArray a
  return r

withIMXDataList :: [IMXData] -> ([MAnyArray] -> IO a) -> IO a
withIMXDataList = mapWith withIMXData

takeIMXData :: MXArray a -> IO IMXData
takeIMXData a = do
  d <- imxData a
  freeMXArray a
  return d

imxFun :: MFun -> IMXFun
imxFun fun a no =
  mapWith withIMXData a $ \m ->
  fun m no >>= mapM takeIMXData

iMXFun :: IMXFun -> MFun
iMXFun fun a no = do
  ia <- mapM imxData a
  fun ia no >>= mapM iMXData

showsApp :: String -> ShowS -> ShowS
showsApp f a s = f ++ '(' : a (')' : s)

{-
showsReshape :: (MIndex,MIndex) -> ShowS -> ShowS
showsReshape (MSubs [],MSubs []) a s = a s
showsReshape (MSubs [],MSubs []) a s = a s
showsReshape r a s = "reshape(" ++ a (',' : shows (realMSize $ mRangeSize r) (')' : s))

showsIMXArrayDelimWith :: (Char,Char,Char) -> (a -> ShowS) -> IMXArray a -> ShowS
showsIMXArrayDelimWith d f a = showsReshape (bounds a) $ showsMListWith d f (elems a)
showsIMXArrayWith :: (a -> ShowS) -> IMXArray a -> ShowS
showsIMXArrayWith = showsIMXArrayDelimWith ('[',',',']')
showsIMXArray :: Show a => IMXArray a -> ShowS
showsIMXArray = showsIMXArrayWith shows
-}

showsMList :: (Char,Char,Char) -> (a -> ShowS) -> [a] -> ShowS
showsMList (l,_,r) _ [] s = l : r : s
showsMList (l,d,r) f (x:xs) s = l : f x (shml xs) where
  shml [] = r : s
  shml (x:xs) = d : f x (shml xs)

showsReshape :: MSize -> (Char -> ShowS) -> ShowS
showsReshape [] f s = f ' ' s
showsReshape [_] f s = f ';' s
showsReshape [1,_] f s = f ',' s
showsReshape z f s = "reshape(" ++ f ' ' (',' : showsMList ('[',' ',']') shows (realMSize z) (')' : s))

showsIMXGenArray :: (Char,Char) -> (a -> ShowS) -> IMXArray a -> ShowS
showsIMXGenArray (l,r) f a = showsReshape (mRangeSize $ bounds a) $ \d -> showsMList (l,d,r) f (elems a)
showsIMXArrayWith :: (a -> ShowS) -> IMXArray a -> ShowS
showsIMXArrayWith = showsIMXGenArray ('[',']')
showsIMXArray :: Show a => IMXArray a -> ShowS
showsIMXArray = showsIMXArrayWith shows

showsComplex :: (Show a, RealFloat a) => Complex a -> ShowS
showsComplex (x:+y) s = "complex(" ++ shows x (',' : shows y (')' : s))

showsMString :: String -> ShowS
showsMString s = showChar '\'' . showString s . showChar '\''

showsIMX :: IMXData -> ShowS
showsIMX IMXNull = showString "[]"
showsIMX (IMXCell a) = showsIMXGenArray ('{','}') shows a
showsIMX (IMXStruct f a) = showsReshape (mRangeSize (r0,r1)) $ \d ->
  showString "struct" . showsMList ('(',',',')') (shf d) (zip f v) where
  shf d (f,v) s = showsMString f $ ',' : showsMList ('{',d,'}') showsIMX v s
  v = transpose $ segment (length f) $ elems a
  ((r0,_),(r1,_)) = bounds a
showsIMX (IMXLogical a) = showsApp "logical"  $ showsIMXArray a
showsIMX (IMXChar a)    = showsApp "char"     $ showsIMXArray a
showsIMX (IMXDouble a)  =                       showsIMXArray a
showsIMX (IMXSingle a)  = showsApp "single"   $ showsIMXArray a
showsIMX (IMXInt8 a)    = showsApp "int8"     $ showsIMXArray a
showsIMX (IMXUint8 a)   = showsApp "uint8"    $ showsIMXArray a
showsIMX (IMXInt16 a)   = showsApp "int16"    $ showsIMXArray a
showsIMX (IMXUint16 a)  = showsApp "uint16"   $ showsIMXArray a
showsIMX (IMXInt32 a)   = showsApp "int32"    $ showsIMXArray a
showsIMX (IMXUint32 a)  = showsApp "uint32"   $ showsIMXArray a
showsIMX (IMXInt64 a)   = showsApp "int64"    $ showsIMXArray a
showsIMX (IMXUint64 a)  = showsApp "uint64"   $ showsIMXArray a
showsIMX (IMXComplexDouble a)   = showsIMXArrayWith showsComplex a
showsIMX (IMXComplexSingle a)   = showsApp "single" $ showsIMXArrayWith showsComplex a
showsIMX (IMXObject c a) = showsApp "class" $ showsIMX a . showChar ',' . showsMString c

instance Show IMXData where
  showsPrec _ = showsIMX

ext: hs

max_stars_repo_name: CornellCAC/haskell-matlab

max_issues_repo_name: CornellCAC/haskell-matlab

max_forks_repo_name: CornellCAC/haskell-matlab


text:
{-|
Module      : HABQTlib.MeasurementProcessing

Functions that deal with optimising and simulating measurements that take place
during tomography.
-}
module HABQTlib.MeasurementProcessing
  ( PurePOVM
  , SingleQbParam
  , measurementProbs
  , svToAngles
  , blochAnglesToSV
  , mkAntipodalPOVM
  , productPOVM
  , simulateMeasuremet
  , optimiseSingleQbPOVM
  ) where

import Control.Applicative (liftA2)
import Control.Newtype.Generics (unpack)
import Data.Complex (mkPolar, polar)
import Data.List (unfoldr)
import Data.Maybe (fromJust)
import qualified Data.Vector as V
import HABQTlib.Data
import HABQTlib.Data.Particle
import qualified Numeric.GSL as GSL
import Numeric.LinearAlgebra (Complex((:+)))
import qualified Numeric.LinearAlgebra as LA
import qualified System.Random.MWC as MWC
import System.Random.MWC.Distributions (categorical)

-- | A POVM consisting of projections onto pure states.
type PurePOVM = V.Vector WeighedPureStateVector

-- | Spherical coordinates of a pure single qubit state on Bloch sphere.
type SingleQbParam = (Double, Double)

-- | Probabilities to measure elements of POVM when performing the measurement
-- over a mixed state.
measurementProbs :: PurePOVM -> DensityMatrix -> V.Vector Double
measurementProbs povm dm = V.map (weighedProb dm) povm

weighedProb :: DensityMatrix -> WeighedPureStateVector -> Double
weighedProb dm (WeighedPureStateVector (w, sv)) = w * pureStateLikelihood sv dm

entropy :: V.Vector Double -> Double
entropy = negate . V.sum . V.map (\p -> p * log p)

povmPointEntropy :: PurePOVM -> DensityMatrix -> Double
povmPointEntropy povm dm = entropy (measurementProbs povm dm)

povmMeanEntropy :: PurePOVM -> ParticleHierarchy -> Double
povmMeanEntropy povm ph =
  let meanEnt = foldOverPts (povmPointEntropy povm) (*) (+) 0
      totalWeight = V.sum . V.map ptsWeight $ ph
      rawEntropy = V.sum . V.map (liftA2 (*) ptsWeight meanEnt) $ ph
   in rawEntropy / totalWeight

-- | Recover a state vector from spherical coordinates.
blochAnglesToSV :: SingleQbParam -> PureStateVector
blochAnglesToSV (th, phi) =
  let z = LA.fromList [1, 0]
      o = LA.fromList [0, 1]
   in PureStateVector . LA.fromColumns . pure $
      LA.scalar (cos (th / 2) :+ 0) * z +
      LA.scalar (mkPolar (sin (th / 2)) phi) * o

-- | Return spherical coordinates of a single qubit pure state (on Bloch
-- sphere).
svToAngles :: PureStateVector -> Maybe SingleQbParam
svToAngles (PureStateVector sv) =
  let s = LA.size sv
      (m0, ph0) = polar $ LA.atIndex sv (0, 0)
      (_, ph1) = polar $ LA.atIndex sv (1, 0)
      ph = ph1 - ph0
      th = 2 * acos m0
   in if s == (2, 1)
        then Just (th, ph)
        else Nothing

-- | Given spherical coordinates, construct a POVM from the given vector and
-- one orthogonal to it.
mkAntipodalPOVM :: SingleQbParam -> PurePOVM
mkAntipodalPOVM c@(th, phi) =
  let sv = blochAnglesToSV c
      sv' = blochAnglesToSV (pi - th, phi + pi)
   in V.fromList $ WeighedPureStateVector <$> [(1, sv), (1, sv')]

reshapeList :: Int -> [a] -> [[a]]
reshapeList n =
  unfoldr
    (\b ->
       if length b < n
         then Nothing
         else Just (splitAt n b))

listToPairs :: [a] -> [(a, a)]
listToPairs = fmap (\(a:[b]) -> (a, b)) . reshapeList 2

pairToList :: (a, a) -> [a]
pairToList (a, b) = [a, b]

-- | Given a list of POVM measurements on sub-systems, construct a POVM over
-- the composite system that includes all of them.
productPOVM :: [PurePOVM] -> PurePOVM
productPOVM sqbPovms =
  let pr ::
           WeighedPureStateVector
        -> WeighedPureStateVector
        -> WeighedPureStateVector
      pr wsv0 wsv1 =
        WeighedPureStateVector (w0 * w1, PureStateVector $ LA.kronecker sv0 sv1)
        where
          up = fmap unpack . unpack
          (w0, sv0) = up wsv0
          (w1, sv1) = up wsv1
   in V.foldl1' (liftA2 pr) $ V.fromList sqbPovms

-- | Approximate most informative separable POVM over a composite system of
-- quantum bits, given a list of single-qubit starting points and a particle
-- distribution.
optimiseSingleQbPOVM ::
     OptIter -- ^ Number of optimisation steps to perform
  -> [PureStateVector] -- ^ single qubit initial states
  -> ParticleHierarchy
  -> PurePOVM
optimiseSingleQbPOVM iter sv0s ph =
  let nq = length sv0s
      method = GSL.NMSimplex2
      precision = 1e-6
      iterations = iter
      initialBox = replicate (2 * nq) (pi / 2)
      estimateDM = getMixedEstimate ph
      obj params =
        negate $ povmPointEntropy povm estimateDM - povmMeanEntropy povm ph
        where
          povm = productPOVM . fmap mkAntipodalPOVM . listToPairs $ params
      start = concatMap (pairToList . fromJust . svToAngles) sv0s
      result = GSL.minimize method precision iterations initialBox obj start
   in productPOVM . fmap mkAntipodalPOVM . listToPairs . fst $ result

-- | Simulate a POVM over a mixed state and return the state vector on which
-- the projection was obtained.
simulateMeasuremet ::
     DensityMatrix -> PurePOVM -> MWC.GenIO -> IO PureStateVector
simulateMeasuremet dm povm gen = do
  let probs = measurementProbs povm dm
      svs = V.map (\(WeighedPureStateVector (_, sv)) -> sv) povm
  idx <- categorical probs gen
  return $ svs V.! idx

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
{-# LANGUAGE ConstraintKinds       #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE LambdaCase            #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE MultiWayIf            #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}

{-|
Stability: alpha

Note that this is an __internal__ module. Use it at your own risk.
-}
module Isumi.Math.FFT.Internal.Polynomial
  (
  -- * Public APIs
    Polynomial(..)
  , polyToVecRep
  , polyFromVecRep
  , polyFromExpCoefRep
  , polyDegree
  , mapPoly
  , UnboxFloatingAppxEq
  , modPoly
  , RealOrComplexPolynomial(..)
  , modPolyG
  , BruunDivisor(..)
  , factorDivisor
  , divToPoly
  , divDegree
  -- * Internal APIs
  , polyMultiply'
  ) where

import           Control.Exception
import           Control.Monad.Loops
import           Control.Monad.ST
import           Data.Complex
import           Data.Foldable                 (for_)
import           Data.STRef
import qualified Data.Vector.Unboxed           as UV
import qualified Data.Vector.Unboxed.Mutable   as MUV
import           Isumi.Math.AppxEq
import           Isumi.Math.FFT.Internal.Utils (findM, isSorted, realToComplex)

-- | A finite polynomial, represented in either vector or list of (exponent,
-- coefficient). Using these constructors implys no check for pre-conditions.
data Polynomial a = PolyVecRep (UV.Vector a)
                  | PolyExpCoefRep [(Int, a)]
                    deriving Show

instance (UV.Unbox a, AppxEq a) => AppxEq (Polynomial a) where
  PolyVecRep v1 ~= PolyVecRep v2 = v1 ~= v2
  PolyExpCoefRep xs ~= PolyExpCoefRep ys = xs ~= ys
  _ ~= _ = False

-- | Construct a polynomial from a vector of coefficients.
--
-- Coefficients of higher degree entries are stored first.
polyFromVecRep :: (UV.Unbox a, Floating a, Eq a)
               => UV.Vector a -> Maybe (Polynomial a)
polyFromVecRep cs = Just . PolyVecRep $ UV.dropWhile (== 0) cs

-- | Construct a polynomial from a list of (exponent, coefficient) pairs
--
-- Higher degree entries are stored first. Unlike vector representation, this
-- representation doesn't need to store every entry, only those have a non-zero
-- coefficient.
polyFromExpCoefRep :: [(Int, a)] -> Maybe (Polynomial a)
polyFromExpCoefRep dcs = if isSorted False (fmap fst dcs)
                            then Just $ PolyExpCoefRep dcs
                            else Nothing

polyToVecRep :: Polynomial a -> Maybe (UV.Vector a)
polyToVecRep (PolyVecRep xs) = Just xs
polyToVecRep _               = Nothing

mapPoly :: (UV.Unbox a, UV.Unbox b)
        => (a -> b)
        -> Polynomial a
        -> Polynomial b
mapPoly f (PolyVecRep v) = PolyVecRep (UV.map f v)
mapPoly f (PolyExpCoefRep ecs) = PolyExpCoefRep $ fmap (\(e, c) -> (e, f c)) ecs

{-| Divisor in the Bruun's FFT algorithm.
'a' is some floating type, such as 'Double'

There are three types of divisors, each of them is represented by a constructor.
-}
data BruunDivisor a =
                      -- | \( x^n - 1 \)
                      BruunDivisorMinus Int
                      -- | \( x^n + cx^{\frac{n}{2}} + 1 \)
                    | BruunDivisorPlus Int a
                      -- | \( x + c \)
                    | BruunDivisorComplex (Complex a)
                      deriving Show

instance (RealFloat a, AppxEq a) => AppxEq (BruunDivisor a) where
  (BruunDivisorMinus n1) ~= (BruunDivisorMinus n2) = n1 ~= n2
  (BruunDivisorPlus n1 c1) ~= (BruunDivisorPlus n2 c2) = n1 ~= n2 && c1 ~= c2
  (BruunDivisorComplex c1) ~= (BruunDivisorComplex c2) = c1 ~= c2
  (BruunDivisorComplex c1) ~= (BruunDivisorMinus n) = n == 1 && c1 ~= (-1)
  x@(BruunDivisorMinus _) ~= y@(BruunDivisorComplex _) = y ~= x
  _ ~= _ = False

divDegree :: BruunDivisor a -> Int
divDegree (BruunDivisorComplex _) = 1
divDegree (BruunDivisorPlus n _)  = n
divDegree (BruunDivisorMinus n)   = n

data RealOrComplexPolynomial a = RealPolynomial (Polynomial a)
                               | ComplexPolynomial (Polynomial (Complex a))

instance (UV.Unbox a, Show a) => Show (RealOrComplexPolynomial a) where
  show (RealPolynomial p)    = show p
  show (ComplexPolynomial p) = show p

{-|
  Convert 'BruunDivisor' to 'Polynomial'
-}
divToPoly :: forall a . (UV.Unbox a, RealFloat a)
          => BruunDivisor a
          -> RealOrComplexPolynomial a
divToPoly = \case
  BruunDivisorMinus deg -> RealPolynomial $
    PolyExpCoefRep [(deg, 1), (0, -1)]
  BruunDivisorPlus deg c -> RealPolynomial $
    PolyExpCoefRep [(deg, 1), (deg `div` 2, c), (0, 1)]
  BruunDivisorComplex c -> ComplexPolynomial $
    PolyExpCoefRep [(1, 1), (0, c)]

type UnboxFloatingAppxEq a = (UV.Unbox a, Floating a, AppxEq a)

modPolyG :: (UnboxFloatingAppxEq a, RealFloat a)
         => RealOrComplexPolynomial a
         -> RealOrComplexPolynomial a
         -> Maybe (RealOrComplexPolynomial a)
modPolyG (RealPolynomial p) (RealPolynomial d) = RealPolynomial <$> modPoly p d
modPolyG (ComplexPolynomial p) (ComplexPolynomial d) =
  ComplexPolynomial <$> modPoly p d
modPolyG (RealPolynomial p) (ComplexPolynomial d) = ComplexPolynomial <$>
  mapPoly realToComplex p `modPoly` d
modPolyG (ComplexPolynomial p) (RealPolynomial d) = ComplexPolynomial <$>
  p `modPoly` mapPoly realToComplex d

{-|
Compute polynomial modulo. Only supports the situation where the first
polynomial is represented by vector,and the second polynomial is represented by
degree-coefficient
-}
modPoly :: UnboxFloatingAppxEq a
        => Polynomial a
        -> Polynomial a
        -> Maybe (Polynomial a)
modPoly _ (PolyExpCoefRep [])               = Nothing
modPoly (PolyVecRep cs) (PolyExpCoefRep ic) = Just . PolyVecRep $ modPoly' cs ic
modPoly _ _                                 = Nothing

{-|
Basically follows the pseudo code from Wikipedia:
  > function n / d:
  >   require d ≠ 0
  >   q ← 0
  >   r ← n       # At each step n = d × q + r
  >   while r ≠ 0 AND degree(r) ≥ degree(d):
  >      t ← lead(r)/lead(d)     # Divide the leading terms
  >      q ← q + t
  >      r ← r − t * d
  >   return (q, r)
-}
modPoly' :: forall a. UnboxFloatingAppxEq a
         => UV.Vector a
         -> [(Int, a)]
         -> UV.Vector a
modPoly' n d = UV.create $ do
  let dDeg = polyDegree . PolyExpCoefRep $ d
  let dHeadCoef :: a = snd . head $ d
  r <- UV.thaw n
  rDeg <- newSTRef (polyDegree . PolyVecRep $ n)
  rIsZero <- ((==0) <$> readSTRef rDeg) >>= newSTRef

  whileM_ ((&&) <$> (not <$> readSTRef rIsZero)
                <*> ((>=dDeg) <$> readSTRef rDeg))
          $ do
    rc <- readSTRef rDeg >>= readCoefAtDeg r
    let tc = rc / dHeadCoef
    tDeg <- (\rDeg' -> rDeg' - dDeg) <$> readSTRef rDeg
    polySubtract' (r, rDeg, rIsZero) (polyMultiply' (tDeg, tc) d)

  rDeg_ <- readSTRef rDeg
  pure $ MUV.drop (MUV.length r - 1 - rDeg_) r

polyDegree :: UV.Unbox a
           => Polynomial a
           -> Int
polyDegree (PolyVecRep v) =
  let len = UV.length v
   in if len == 0
         then 0
         else len - 1
polyDegree (PolyExpCoefRep ic) =
  case ic of
    []    -> 0
    (x:_) -> fst x

readCoefAtDeg :: UV.Unbox a => UV.MVector s a -> Int -> ST s a
readCoefAtDeg v d = MUV.read v (MUV.length v - 1 - d)

updateCoefAtDeg :: UV.Unbox a => UV.MVector s a -> Int -> a -> ST s ()
updateCoefAtDeg v d = MUV.write v (MUV.length v - 1 - d)

polyMultiply' :: Num a => (Int, a) -> [(Int, a)] -> [(Int, a)]
polyMultiply' (d, c) = fmap (\(d', c') -> (d + d', c * c'))

polySubtract' :: UnboxFloatingAppxEq a
              => (UV.MVector s a, STRef s Int, STRef s Bool)
              -> [(Int, a)]
              -> ST s ()
polySubtract' (r, rDeg, rIsZero) dcs = do
  for_ dcs $ \(d, c) -> do
    c' <- (\x -> x - c) <$> readCoefAtDeg r d
    updateCoefAtDeg r d c'
  -- update 'rDeg' and 'rIsZero'
  -- search from rDeg to lower degrees, until a non-zero coefficient is found.
  rDeg_ <- readSTRef rDeg
  degM <- findM (fmap (~/= 0) . readCoefAtDeg r)
            [rDeg_, (rDeg_-1) .. 0]
  case degM of
    Nothing -> do
      writeSTRef rDeg 0
      writeSTRef rIsZero True
    Just d -> writeSTRef rDeg d

-- | Factor a divisor into two divisors.
--
-- If current divisor is of degree 1, it can not be factored again,
-- then 'Nothing' will be returned.
factorDivisor :: RealFloat a
              => BruunDivisor a
              -> Maybe (BruunDivisor a, BruunDivisor a)
factorDivisor (BruunDivisorComplex _) = Nothing
factorDivisor (BruunDivisorMinus n) =
  if n == 1
     then Nothing
     else
       let n' = n `div` 2
        in Just (BruunDivisorMinus n', BruunDivisorPlus n' 0)
factorDivisor (BruunDivisorPlus n c) =
  if | n == 1 -> Nothing
     | n `rem` 4 == 0 -> Just $ factorDPToDPs n c
     | otherwise -> Just $ factorDPToDCs c

-- | Factor divisor in plus form into two divisors in plus form.
-- Doesn't check for pre-condition.
factorDPToDPs :: Floating a => Int -> a -> (BruunDivisor a, BruunDivisor a)
factorDPToDPs n c =
  let n' = n `div` 2
      c' = sqrt (2 - c)
   in assert (n' >= 2)
        (BruunDivisorPlus n' c', BruunDivisorPlus n' (-c'))

-- | Factor divisor in plus form into two divisors in complex form.
-- Doesn't check for pre-condition.
factorDPToDCs :: RealFloat a
              => a
              -> (BruunDivisor a, BruunDivisor a)
factorDPToDCs c =
  (BruunDivisorComplex (c' / 2 + d), BruunDivisorComplex (c' / 2 - d))
  where
    c' = realToComplex c
    d = sqrt (c' * c' - 4) / 2

max_stars_repo_name: IsumiF/fft-bruun

max_issues_repo_name: IsumiF/fft-bruun

max_forks_repo_name: IsumiF/fft-bruun


text:
{-|

 Module      :  HF.Gauss
 Copyright   :  Copyright (c) David Schlegel
 License     :  BSD
 Maintainer  :  David Schlegel
 Stability   :  experimental
 Portability :  Haskell

Gaussian Integral evaluation plays an important role in quantum chemistry. Here functions will be provided to compute the most important integrals involving gaussian-type orbitals (GTOs).

An unnormalized primitive cartesian Gaussian Orbital has the form

<<centered_gaussian.svg centered_gaussian>>

where /A/ denotes the center of the orbital.

An unnormalized contracted Gaussian Orbital is a linear combination of primitive Gaussian:

<<centered_contracted_gaussian.svg centered_contracted_gaussian>>

"Gauss" uses the datatype contstructors 'PG' (Primitive Gaussian) and 'Ctr' (Contraction), defined in "HF.Data".

Notice: Formulas of Gaussian Integrals are taken from:
	Fundamentals of Molecular Integral Evaluation
	by Justin T. Fermann and Edward F. Valeev

-}

module HF.Gauss (
-- * Basic functions
distance, frac, center, itersum, erf, f, factorial2, general_erf,
-- * Normalization
normCtr, normPG,
-- * Integral Evaluation
s_12, t_12,
-- * Contraction operations
zipContractionWith, constr_matrix

) where

import HF.Data
import Numeric.Container hiding (linspace)
import Numeric.LinearAlgebra.Data hiding (linspace)
import Data.Maybe
import Control.DeepSeq
import Debug.Trace

---------------------
---------------------
--Helper Functions---
---------------------
---------------------

itersum function range
	| range == [] = 0
	| otherwise = (function $ head range) + itersum function (tail range)

-- |Calculates the factorial f(x) = x!
factorial n
    | n < 0     = 0
    | otherwise = product [1..n]

-- |Calculates the Doublefactorial f(x) = x!!
factorial2 :: (Show a, Eq a, Num a) => a -> a
factorial2 0 = 1
factorial2 1 = 1
factorial2 (-1) = 1
factorial2 n = n * factorial2 (n-2)

-- |Calculates the binomialcoefficient n over k
--binom :: Integer -> Integer -> Double
binom n 0 = 1
binom 0 k = 0
binom n k = binom (n-1) (k-1) * n `div` k

-- |Computes the distance between to vectors:
--
-- <<norm.svg norm>>
distance :: Vector Double -> Vector Double -> Double
distance rA rB = norm2(rA - rB)

-- | Calculates effective coefficient:
--
-- >>> frac alpha beta = alpha*beta/(alpha + beta)
frac ::  Double -> Double -> Double
frac alpha beta = alpha*beta/(alpha + beta)

-- | Calculates weighted radius of two gaussians
center :: Double -- ^ alpha
			-> Vector Double -- ^ rA
			-> Double -- ^ beta
			-> Vector Double -- ^ rB
			-> Vector Double
center a rA b rB = ((scalar b * rB) * (scalar a * rA))/scalar (a+b)

-- |error function calculation
erf :: Double -> Double
erf x = 2/sqrt pi * sum [x/(2*n+1) * product [-x^2/k |  k <- [1..n] ] | n <- [0..100]]

--error function calculation
f_0 :: Double -> Double
f_0 = erf

-- | General error function. See Eq. 4.118 Thijssen.
general_erf nu u
	| nu == 0 = sqrt pi * erf (sqrt u) * 1/(2*sqrt u)
	| otherwise = 1/(2*u) * (2*nu - 1) * (general_erf (nu-1) u ) - exp (-u)

--------------------------
---Important Functions----
--------------------------

-- |Calculates normalization factor of contracted Gaussian with arbitrary angular momentum
normCtr :: Ctr -> Double
normCtr contr = (prefactor * summation)**(-1.0/2.0)
	--See also Eq. 2.11
	where
		(l, m, n) = lmncontract contr
		--(l,m,n) = (fromIntegral l_, fromIntegral m_, fromIntegral n_) --This looks quite dirty
		n_sum = lengthcontract contr -1
		a = toList $ coefflist contr
		alp = [alpha prim | prim <- (gaussians contr)]
		mom = fromIntegral $ momentumctr contr
		prefactor = 1.0/(2.0**mom)*pi**(3.0/2.0)* factorial2 (2*l-1) * factorial2 (2*m-1)* factorial2 (2*n-1)
		summation = sum $ concat $ [[(a !! i)*(a !! j)/((alp !! i +alp !! j)**(mom + 3.0/2.0)) | i <-[0..n_sum]]| j <- [0..n_sum]]

-- |Calculates normalization factor of primitive Gaussian with arbitrary angular momentum
normPG :: PG -> Double
normPG (PG lmn alpha pos) = (2*alpha/pi)**(3/4) * (a/b)**(1/2)
	where
		(l, m, n) = (fromIntegral $ lmn !! 0, fromIntegral $ lmn !! 1, fromIntegral $ lmn !! 2)
		a = (8*alpha)**(l+m+n) * factorial l * factorial m * factorial n
		b = factorial (2*l) * factorial (2*m) * factorial (2*n)

-- |Calculates f_k (l1 l2 pa_x pb_x) used in Gaussian Product Theorem
f :: Int -> Int -> Int -> Double -> Double -> Double
--trace ("f" ++ " " ++  show k ++ " " ++ show l1_ ++ " " ++ show l2_ ++ " " ++ show pa_x ++ " " ++ show pb_x) $
f k l1_ l2_ pa_x pb_x = sum $ concat $ [[pa_x**(fromIntegral (l1_-i))  *  pb_x**(fromIntegral(l2_-j)) * fromIntegral ((binom l1_ i) * (binom l2_ j))| i <- [0..l1_], i+j <= k]| j <-[0..l2_]]
--See also Eq. 2.45

-- |Evaluates a given function for two contractions. The operation is symmetric in contraction arguments.
zipContractionWith :: (PG -> PG -> Double) -- ^ Function of two primitive Gaussians
							-> Ctr -- ^ Contraction
							-> Ctr -- ^ Contraction
							-> Double
zipContractionWith zipfunction (Ctr pglist1 coeffs1) (Ctr pglist2 coeffs2) = pr1 * pr2 * value
	where
		coefflist1 = toList coeffs1
		coefflist2 = toList coeffs2
		n1 = length coefflist1 -1
		n2 = length coefflist2 -1
		pr1 = normCtr (Ctr pglist1 coeffs1) --I am not so sure about this
		pr2 = normCtr (Ctr pglist2 coeffs2) --I am not so sure about this
		value = sum $ [(coefflist1 !! i)* (coefflist2 !! j)* (zipfunction (pglist1 !! i) (pglist2 !! j)) | i <- [0..n1], j <- [0..n2]]

-- |Construct a matrix out of a list of Contractions and a function for two primitive gaussians
constr_matrix :: [Ctr] -- ^ List of Contraction
						-> (PG -> PG -> Double) -- ^ Function of two primitive Gaussians
						-> Matrix Double -- ^ Matrix of dimensions (length [Ctr]) x (length [Ctr])
constr_matrix contractionlist function = buildMatrix (length contractionlist) (length contractionlist) (\(i,j) -> zipContractionWith function (contractionlist !! i) (contractionlist !! j) )

-- |Calculates overlap integral of two primitive gaussians

-- | <<s_12.svg s_12>>
s_12 :: PG -> PG -> Double
s_12 (PG lmn1 alpha1 pos1) (PG lmn2 alpha2 pos2) = trace ("s_12:   range " ++ show range  ++ "   f 0 =" ++ show (factor 0) ++ "   f 1 = " ++ show (factor 1)++ "   f 2 =" ++ show (factor 2)) $ prefactor * (factor 0) * (factor 1) * (factor 2) --pref * I_x * I_y * I_z
	where
		(l1, m1, n1) = (lmn1 !! 0, lmn1 !! 1, lmn1 !! 2)
		(l2, m2, n2) = (lmn2 !! 0, lmn2 !! 1, lmn2 !! 2)
		g = alpha1 + alpha2
		prefactor = exp (-alpha1*alpha2 * (distance pos1 pos2) /g)
		p = center alpha1 pos1 alpha2 pos2
		pa = toList $ p  - pos1
		pb = toList $ p  - pos2
		range = [0..(fromIntegral (l1+l2) /2)] :: (Enum a, Fractional a) => [a]
		-- See also Eq. 3.15
		function k i = force $ f (round (2*i)) l1 l2 (pa !! k) (pb !! k) *  (factorial2 (2.0* i -1.0)) * ((2*g)** (- i)) * (pi/g)**(1/2)
		factor k = itersum (function k) range
		--i k = sum $ [f (round (2*i)) l1 l2 (pa !! k) (pb !! k) * factorial2 (2.0* i -1.0) * ((2*g)** (- i)) * (pi/g)**(1/2) | i <- range ]

-- |Calculates kinetic energy integral of two primitive gaussians

-- | <<t_12.svg t_12>>
t_12 :: PG -> PG -> Double
t_12 (PG lmn1 alpha1 pos1) (PG lmn2 alpha2 pos2) =  (i 0) + (i 1) + (i 2) -- I_x + I_y + I_z
	where
		--lifts or lowers l, m or n, indicated by index k for a lmn-list, respectively
		pp1 list k = [if i == k then list !! i + 1 else list !! i | i<-[0..2]]	--lifts
		mm1 list k = [if i == k then list !! i - 1 else list !! i | i<-[0..2]]	--lowers
		--overlaps where parts of l, m or n ar lifted +1 or lowered -1
		m1m1 k = s_12 (PG (mm1 lmn1 k) alpha1 pos1) (PG (mm1 lmn2 k) alpha2 pos2) -- <-1|-1>_k
		p1p1 k = s_12 (PG (pp1 lmn1 k) alpha1 pos1) (PG (pp1 lmn2 k) alpha2 pos2) -- <+1|+1>_k
		p1m1 k = s_12 (PG (pp1 lmn1 k) alpha1 pos1) (PG (mm1 lmn2 k) alpha2 pos2) -- <+1|-1>_k
		m1p1 k = s_12 (PG (mm1 lmn1 k) alpha1 pos1) (PG (pp1 lmn2 k) alpha2 pos2) -- <-1|+1>_k
		--See also Eq.4.13
		i k = 1/2 * (fromIntegral (lmn1 !! k)) * (fromIntegral (lmn2 !! k)) * (m1m1 k)
			+ 2 * alpha1 * alpha2 * (p1p1 k)
			- alpha1 * (fromIntegral (lmn2 !! k)) * (p1m1 k)
			- alpha2 * (fromIntegral (lmn1 !! k)) * (m1p1 k)

{-
-- |Calculates nuclear attraction integral of two primitive gaussians
v_12 :: Double -> Vecotr Double -> PG -> PG -> Double
v_12 z r_c (PG lmn1 alpha1 pos1) (PG lmn2 alpha2 pos2) = z*pi*pre/g * 2* summation
	where
		ax, ay, az = (fromList pos1) !! 0, (fromList pos1) !! 1, (fromList pos1) !! 2
		bx, by, bz = (fromList pos2) !! 0, (fromList pos2) !! 1, (fromList pos2) !! 2
		r_cx, r_cy, r_cz = (fromList r_c) !! 0, (fromList r_c) !! 1, (fromList r_c) !! 2
		pref = exp (-alpha1*alpha2 * (distance pos1 pos2) /g)
		g = alpha1 + alpha2
		p = 1/g * (alpha1*pos1 + alpha2*pos2)
		px, py, pz = (fromList p) !! 0, (fromList p) !! 1, (fromList p) !! 2
		eta = (alpha1 + alpha2)/g
		mu_x = l1 + l2 - 2*(ijk1 + ijk2) - (opq1 + opq2)
		mu_y = l1 + l2 - 2*(j1 + j2) - (p1 + p2)
		mu_z = l1 + l2 - 2*(k1 + k2) - (q1 + q2)
		nu = mu_x + mu_y + mu_z - (u + v + w)
		summation = general_erf nu (g * distance p r_c) * a_x * a_y * a_z

		--ranges:
		ijk1_range = [0..floor $ (fromIntegral l1 /2)] :: (Integral a, Enum a) => [a]
		ijk2_range = [0..floor $ (fromIntegral l2 /2)] :: (Integral a, Enum a) => [a]
		opq1_range ijk1 = [0..floor $ (fromIntegral l1 - 2*ijk1)] :: (Integral a, Enum a) => [a]
		opq2_range ijk2 = [0..floor $ (fromIntegral l2 - 2*ijk2)] :: (Integral a, Enum a) => [a]
		rst_range opq1 opq2 = [0..floor $ fromIntegral (opq1+l2) /2] :: (Integral a, Enum a) => [a]
		u_range mu_x = [0..floor $ (fromIntegral mu_x /2)] :: (Integral a, Enum a) => [a]

		a_x = (-1)**(l1+l2) * (factorial l1) * (factorial l2) *
		(-1)**(opq2 + r) * factorial (opq1 + opq2) * 1/(4**(ijk1+ijk2+r) * factorial ijk1 * factorial ijk2 * factorial opq1 * factorial opq2 * factorial r) *
		alpha1**(opq2-ijk1-r) * alpha2**(opq1-ijk2-r) * (ax-bx)**(opq1+opq2-2*r) * 1/(factorial (l1-2*ijk1-opq1) * factorial (l2-2*ijk2-opq2) * factorial (opq1-2*opq2-2*r)) *
		* [ (-1)**u * factorial mu_x * (px -r_cx)**(mu_x - 2*u) * 1/ ((4)**u * factorial u * factorial (mu_x-2*u) * g**(opq1 + opq2 - r + u)) | u <- [0..1/2 * mu_x]]
-}

--v_12_test ::  Double -> Vector Double -> PG -> PG -> [Double]
{- v_12_test z r_c (PG lmn1 alpha1 pos1) (PG lmn2 alpha2 pos2) = summation 0
	where
		--(ax, ay, az) = ((toList pos1) !! 0, (toList pos1) !! 1, (toList pos1) !! 2)
		--(bx, by, bz) = ((toList pos2) !! 0, (toList pos2) !! 1, (toList pos2) !! 2)
		--(r_cx, r_cy, r_cz) = ((toList r_c) !! 0, (toList r_c) !! 1, (toList r_c) !! 2)
		--lmn1_ = fromList lmn1
		--lmn2_ = fromList lmn2

		pref = exp (-alpha1*alpha2 * (distance pos1 pos2) /g)
		g = alpha1 + alpha2
		p = scale (1/g) (add (scale alpha1 pos1) (scale alpha2 pos2))
		--(px, py, pz) = ((toList p) !! 0, (toList p) !! 1, (toList p) !! 2)

		--f_nu ijk1 ijk2 opq1 opq2 = general_erf (nu ijk1 ijk2 opq1 opq2) (g * distance p r_c)

		eta = (alpha1 + alpha2)/g
		mu ijk1 ijk2 opq1 opq2 s = fromIntegral (atIndex lmn1 s) + fromIntegral (atIndex lmn2 s) - 2*(ijk1 + ijk2) - (opq1 + opq2)
		--nu ijk1 ijk2 opq1 opq2 = (mu ijk1 ijk2 opq1 opq2 s) + (mu ijk1 ijk2 opq1 opq2 s) + (mu ijk1 ijk2 opq1 opq2 s) - (u v q)

		--ranges alle ranges (X,Y,Z) mit einer klappe schlagen, dafür index s:
		ijk1_range s = [0..fromIntegral $ floor $ fromIntegral (lmn1 !! s) /2] :: [Double]
		ijk2_range s = [0..fromIntegral $ floor $ fromIntegral (lmn2 !! s) /2] :: [Double]
		opq1_range ijk1 s = [0..fromIntegral $ floor $ (fromIntegral (lmn1 !! s) - 2*ijk1)] :: [Double]
		opq2_range ijk2 s = [0..fromIntegral $ floor $ (fromIntegral (lmn2 !! s) - 2*ijk2)] :: [Double]
		rst_range opq1 opq2 s = [0..fromIntegral $ floor $ (opq1+opq2) /2] :: [Double]
		uvw_range ijk1 ijk2 opq1 opq2 s = [0..fromIntegral $ floor $ ((mu ijk1 ijk2 opq1 opq2 s) /2)] :: [Double]

		--a_x :: Num a => a -> a -> a -> a -> a -> Double

		--(-1)**(l1+l2) * (factorial fromIntegral (lmn1 !! s)) * (factorial fromIntegral (lmn2 !! s)) *
		function ijk1 ijk2 opq1 opq2 rst s = (-1)**(opq2 + rst) * factorial (opq1 + opq2) * 1/(4**(ijk1+ijk2+rst) * factorial ijk1 * factorial ijk2 * factorial opq1 * factorial opq2 * factorial rst) * alpha1**(opq2-ijk1-rst) * alpha2**(opq1-ijk2-rst) * (atIndex pos1 s - atIndex pos2 s)**(opq1+opq2-2*rst) * 1/(factorial ((lmn1 !! s)-2*ijk1-opq1) * factorial ((lmn2 !! s)-2*ijk2-opq2) * factorial (opq1-2*opq2-2*rst)) * sum $ [ (-1)**uvw * factorial (mu ijk1 ijk2 opq1 opq2 s) * ((atIndex p s) - (atIndex r_c s))**((mu ijk1 ijk2 opq1 opq2 s) - 2*uvw) * 1/ ((4)**uvw * factorial uvw * factorial ((mu ijk1 ijk2 opq1 opq2 s)-2*uvw) * g**(opq1 + opq2 - rst + uvw)) | uvw <- (uvw_range ijk1 ijk2 opq1 opq2 s)]

		summation s = (-1)**(lmn1 !! s + lmn2 !! s) * (factorial fromIntegral (lmn1 !! s)) * (factorial fromIntegral (lmn2 !! s)) * sum $ concat $ concat $ [[[function ijk1 ijk2 opq1 opq2 rst s | rst <- (rst_range opq1 opq2 s)] | opq1 <- (opq1_range ijk1 s), opq2 <- (opq2_range ijk2 s)]| ijk1 <- ijk1_range s, ijk2 <- ijk2_range s]
-}

-----------------
--Deprecated-----
-----------------

-- Here all importand integrals involving gaussian functions will be evaluated.

--Overlap integral for s orbitals
--calculates <1s,alpha, A | 1s, beta, B>
overlaps :: Double -> Double -> Vector Double -> Vector Double -> Double
overlaps alpha beta rA rB = prefactor * exp exponent
	where 	prefactor = (pi/(alpha + beta))**(3/2)
		exponent = - (frac alpha beta * distance rA rB)

--Kinetic integral
--calculates <1s,alpha, A | - Laplace | 1s, beta, B>
kinetic :: Double -> Double -> Vector Double -> Vector Double -> Double
kinetic alpha beta rA rB = prefactor * exp(exponent)
	where 	prefactor = (frac alpha beta)
			 *(6 -4 *(frac alpha beta) * distance rA rB)
			 *(pi/(alpha + beta))**(3/2)
		exponent = - ((frac alpha beta) * distance rA rB)

--Nuclear interaction integral
--calculates <1s,alpha, A | - Z/r_C | 1s, beta, B>
nuclear :: Double -> Double -> Vector Double -> Vector Double -> Vector Double -> Double -> Double
nuclear alpha beta rA rB rC z = pref1 * pref2 * exp(exponent)
	where 	pref1 = -2*pi*z/(alpha + beta)
		pref2 = f_0 arg
		arg   = (distance (center alpha rA beta rB) rC )   * (alpha*beta)
		exponent = - ((frac alpha beta) * (distance rA rB))

--two-electron integral
--important!
--calculates <1s,alpha, A ; 1s, beta, B | 1s,gamma, C ; 1s, delta, D >
twoelectron :: Double -> Double -> Double -> Double -> Vector Double -> Vector Double -> Vector Double -> Vector Double -> Double
twoelectron a b g d rA rB rC rD = pref * exp(exponent) * (f_0 arg)
	where	pref = 			(2*pi**(5/2))/
				( (a+g)*(b+d)*(a+b+g+d)**(1/2)   )
		exponent = 	-a*g/(a+g) * (distance rA rC)
				-b*d/(b+d) * (distance rB rD)
		arg = (a+g)*(b+d)/(a+b+g+d) * ( distance (center a rA g rC) (center b rB d rD) )

--Overlap Matrix S_pq = <p|q>
--overlap p q  = t
--where
--	t = buildMatrix (n*n) (n*n) (\(i,j) -> overlaps alpha beta rA rB )

max_stars_repo_name: davidschlegel/hartree-fock

max_issues_repo_name: davidschlegel/hartree-fock

max_forks_repo_name: davidschlegel/hartree-fock


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

module Utils where

import qualified Graphics.GL as GLRaw
import qualified Graphics.Rendering.OpenGL.GL as GL
import           Numeric.LinearAlgebra ((><), Vector, Matrix, scale, cross, dot, tr)
import qualified Numeric.LinearAlgebra as L
import           Data.Vector.Storable ((!))
import qualified Numeric.LinearAlgebra.Data as LD
import qualified Data.Vector.Storable as VS

step :: Float
step = 0.02

norm :: Vector Float -> Float
norm = realToFrac . L.norm_2

posToMatV :: Vector Float -> Vector Float
posToMatV v =
            [1,0,0,0
             ,0,1,0,0
             ,0,0,1,0
             ,VS.unsafeIndex v 0,VS.unsafeIndex v 1, VS.unsafeIndex v 2,1]

posToMat = tr . LD.reshape 4 . posToMatV

glTexture :: GL.GLuint -> GLRaw.GLenum
glTexture = (GLRaw.GL_TEXTURE0+) . fromIntegral

lookAt' :: Vector Float -> Vector Float -> Vector Float -> Matrix Float
lookAt' up eye at =
  let zaxis = normalize $ eye - at
      xaxis = normalize $ cross zaxis up
      yaxis = cross xaxis zaxis
      crd x = negate $ dot eye x
  in (4><4) [ xaxis!0,   yaxis!0,   zaxis!0,   0
            , xaxis!1,   yaxis!1,   zaxis!1,   0
            , xaxis!2,   yaxis!2,   zaxis!2,   0
            , crd xaxis, crd yaxis, crd zaxis, 1]

pointAt :: Vector Float -> Vector Float
pointAt sp =
  let up =  [0,1,0]
      eye = [0,0,0]
      at  = sp
      zaxis = normalize $ eye - at
      xaxis = normalize $ cross zaxis up
      yaxis = cross xaxis zaxis
      crd x = negate $ dot eye x
  in [ xaxis!0,   yaxis!0,   zaxis!0,   0
     , xaxis!1,   yaxis!1,   zaxis!1,   0
     , xaxis!2,   yaxis!2,   zaxis!2,   0
     , crd xaxis, crd yaxis, crd zaxis, 1]

lookAt :: Vector Float -> Vector Float -> Matrix Float
lookAt = lookAt' [0,1,0]

normalize :: Vector Float -> Vector Float
normalize v =
  let n = (1/) . norm $ v
  in scale n v

remTrans = mat3ToMat4 . mat4ToMat3

mat4ToMat3 :: Matrix Float -> Matrix Float
mat4ToMat3 = LD.takeRows 3 . LD.takeColumns 3

mat3ToMat4 :: Matrix Float -> Matrix Float
mat3ToMat4 = helper . LD.flatten
  -- LD.fromLists . foldr (\a -> ((a ++ [0]):)) [[0,0,0,1]] . LD.toLists
  where
    helper v = (4><4) [ VS.unsafeIndex v 0, VS.unsafeIndex v 1, VS.unsafeIndex v 2, 0
                      , VS.unsafeIndex v 3, VS.unsafeIndex v 4, VS.unsafeIndex v 5, 0
                      , VS.unsafeIndex v 6, VS.unsafeIndex v 7, VS.unsafeIndex v 8, 0
                      ,                  0,                  0,                  0, 1]

projectionMatrix :: Matrix Float
projectionMatrix =
  (4><4) [ t / as, 0,           0, 0
         , 0,      t,           0, 0
         , 0,      0, (f+n)/(n-f),(2*f*n)/(n-f)
         , 0,      0,          -1, 0]
  where
    n = 0.1
    f = 200
    t = 1 / tan ((pi/3) / 2)
    as= 1366/768

maxSpeed = 10

clampSpeed :: Vector Float -> Vector Float
clampSpeed v =
  if norm v >= maxSpeed
  then LD.cmap (*maxSpeed) . normalize $ v
  else v

max_stars_repo_name: Antystenes/CPG

max_issues_repo_name: Antystenes/CPG

max_forks_repo_name: Antystenes/CPG


text:
{-# LANGUAGE BangPatterns , RankNTypes, GADTs, DataKinds #-}

{- | The 'Numerical.HBLAS.BLAS.Level2' module provides a fully general
yet type safe Level2 BLAS API.

When in doubt about the semantics of an operation,
consult your system's BLAS api documentation, or just read the documentation
for
<https://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm the Intel MKL BLAS distribution>

A few basic notes about how to invoke BLAS routines.

Many BLAS operations take one or more arguments of type 'Transpose'.
'Tranpose' has the following different constructors, which tell BLAS
routines what transformation to implicitly apply to an input matrix @mat@ with dimension @n x m@.

*  'NoTranspose' leaves the matrix @mat@ as is.

* 'Transpose' treats the @mat@ as being implicitly transposed, with dimension
    @m x n@. Entry @mat(i,j)@ being treated as actually being the entry
    @mat(j,i)@. For Real matrices this is also the matrix adjoint operation.
    ie @Tranpose(mat)(i,j)=mat(j,i)@

*  'ConjNoTranspose' will implicitly conjugate @mat@, which is a no op for Real ('Float' or 'Double') matrices, but for
'Complex Float' and 'Complex Double' matrices, a given matrix entry @mat(i,j)==x':+'y@
will be treated as actually being  @conjugate(mat)(i,j)=y':+'x@.

* 'ConjTranpose' will implicitly transpose and conjugate the input matrix.
ConjugateTranpose acts as matrix adjoint for both real and complex matrices.

The *gemm operations  work as follows (using 'sgemm' as an example):

* @'sgemm trLeft trRight alpha beta left right result'@, where @trLeft@ and @trRight@
are values of type 'Transpose' that respectively act on the matrices @left@ and @right@.

* the generalized matrix computation thusly formed can be viewed as being
@result = alpha * trLeft(left) * trRight(right) + beta * result@

the *gemv operations are akin to the *gemm operations, but with @right@ and @result@
being vectors rather than matrices.

the *trsv operations solve for @x@ in the equation @A x = y@ given @A@ and @y@.
The 'MatUpLo' argument determines if the matrix should be treated as upper or
lower triangular and 'MatDiag' determines if the triangular solver should treat
the diagonal of the matrix as being all 1's or not.  A general pattern of invocation
would be @'strsv' matuplo  tranposeMatA  matdiag  matrixA  xVector@.
A key detail to note is that the input vector is ALSO the result vector,
ie 'strsv' and friends updates the vector place.

-}

module Numerical.HBLAS.BLAS.Level2(
        sgbmv
        ,dgbmv
        ,cgbmv
        ,zgbmv

        ,sgemv
        ,dgemv
        ,cgemv
        ,zgemv

        ,sger
        ,dger
        ,cgerc
        ,zgerc
        ,cgeru
        ,zgeru

        ,chbmv
        ,zhbmv

        ,chemv
        ,zhemv

        ,cher
        ,zher
        ,cher2
        ,zher2

        ,chpmv
        ,zhpmv

        ,chpr
        ,zhpr
        ,chpr2
        ,zhpr2

        ,ssbmv
        ,dsbmv

        ,sspmv
        ,dspmv

        ,sspr
        ,dspr
        ,sspr2
        ,dspr2

        ,ssymv
        ,dsymv

        ,ssyr
        ,dsyr
        ,ssyr2
        ,dsyr2

        ,stbmv
        ,dtbmv
        ,ctbmv
        ,ztbmv

        ,stbsv
        ,dtbsv
        ,ctbsv
        ,ztbsv

        ,stpmv
        ,dtpmv
        ,ctpmv
        ,ztpmv

        ,stpsv
        ,dtpsv
        ,ctpsv
        ,ztpsv

        ,strmv
        ,dtrmv
        ,ctrmv
        ,ztrmv

        ,strsv
        ,dtrsv
        ,ctrsv
        ,ztrsv
            ) where

import Numerical.HBLAS.UtilsFFI
import Numerical.HBLAS.BLAS.FFI.Level2
import Numerical.HBLAS.BLAS.Internal.Level2
import Control.Monad.Primitive
import Data.Complex

sgbmv :: PrimMonad m => GbmvFun Float orient (PrimState m) m
sgbmv = gbmvAbstraction "sgbmv" cblas_sgbmv_safe cblas_sgbmv_unsafe (\x f -> f x)

dgbmv :: PrimMonad m => GbmvFun Double orient (PrimState m) m
dgbmv = gbmvAbstraction "dgbmv" cblas_dgbmv_safe cblas_dgbmv_unsafe (\x f -> f x)

cgbmv :: PrimMonad m => GbmvFun (Complex Float) orient (PrimState m) m
cgbmv = gbmvAbstraction "cgbmv" cblas_cgbmv_safe cblas_cgbmv_unsafe withRStorable_

zgbmv :: PrimMonad m => GbmvFun (Complex Double) orient (PrimState m) m
zgbmv = gbmvAbstraction "zgbmv" cblas_zgbmv_safe cblas_zgbmv_unsafe withRStorable_

sgemv :: PrimMonad m => GemvFun Float orient (PrimState m) m
sgemv = gemvAbstraction "sgemv" cblas_sgemv_safe cblas_sgemv_unsafe (flip ($))

dgemv :: PrimMonad m => GemvFun Double orient (PrimState m) m
dgemv = gemvAbstraction "dgemv" cblas_dgemv_safe cblas_dgemv_unsafe (flip ($))

cgemv :: PrimMonad m => GemvFun (Complex Float) orient (PrimState m) m
cgemv = gemvAbstraction "cgemv" cblas_cgemv_safe cblas_cgemv_unsafe withRStorable_

zgemv :: PrimMonad m => GemvFun (Complex Double) orient (PrimState m) m
zgemv = gemvAbstraction "zgemv" cblas_zgemv_safe cblas_zgemv_unsafe withRStorable_

sger :: PrimMonad m => GerFun Float orient (PrimState m) m
sger = gerAbstraction "sger" cblas_sger_safe cblas_sger_unsafe (\x f -> f x)

dger :: PrimMonad m => GerFun Double orient (PrimState m) m
dger = gerAbstraction "dger" cblas_dger_safe cblas_dger_unsafe (\x f -> f x)

cgerc :: PrimMonad m => GerFun (Complex Float) orient (PrimState m) m
cgerc = gerAbstraction "cgerc" cblas_cgerc_safe cblas_cgerc_unsafe withRStorable_

zgerc :: PrimMonad m => GerFun (Complex Double) orient (PrimState m) m
zgerc = gerAbstraction "zgerc" cblas_zgerc_safe cblas_zgerc_unsafe withRStorable_

cgeru :: PrimMonad m => GerFun (Complex Float) orient (PrimState m) m
cgeru = gerAbstraction "cgeru" cblas_cgeru_safe cblas_cgeru_unsafe withRStorable_

zgeru :: PrimMonad m => GerFun (Complex Double) orient (PrimState m) m
zgeru = gerAbstraction "zgeru" cblas_zgeru_safe cblas_zgeru_unsafe withRStorable_

chbmv :: PrimMonad m => HbmvFun (Complex Float) orient (PrimState m) m
chbmv = hbmvAbstraction "chbmv" cblas_chbmv_safe cblas_chbmv_unsafe withRStorable_

zhbmv :: PrimMonad m => HbmvFun (Complex Double) orient (PrimState m) m
zhbmv = hbmvAbstraction "zhbmv" cblas_zhbmv_safe cblas_zhbmv_unsafe withRStorable_

chemv :: PrimMonad m => HemvFun (Complex Float) orient (PrimState m) m
chemv = hemvAbstraction "chemv" cblas_chemv_safe cblas_chemv_unsafe withRStorable_

zhemv :: PrimMonad m => HemvFun (Complex Double) orient (PrimState m) m
zhemv = hemvAbstraction "zhemv" cblas_zhemv_safe cblas_zhemv_unsafe withRStorable_

cher :: PrimMonad m => HerFun Float (Complex Float) orient (PrimState m) m
cher = herAbstraction "cher" cblas_cher_safe cblas_cher_unsafe (\x f -> f x)

zher :: PrimMonad m => HerFun Double (Complex Double) orient (PrimState m) m
zher = herAbstraction "zher" cblas_zher_safe cblas_zher_unsafe (\x f -> f x)

cher2 :: PrimMonad m => Her2Fun (Complex Float) orient (PrimState m) m
cher2 = her2Abstraction "cher2" cblas_cher2_safe cblas_cher2_unsafe withRStorable_

zher2 :: PrimMonad m => Her2Fun (Complex Double) orient (PrimState m) m
zher2 = her2Abstraction "zher2" cblas_zher2_safe cblas_zher2_unsafe withRStorable_

chpmv :: PrimMonad m => HpmvFun (Complex Float) orient (PrimState m) m
chpmv = hpmvAbstraction "chpmv" cblas_chpmv_safe cblas_chpmv_unsafe withRStorable_

zhpmv :: PrimMonad m => HpmvFun (Complex Double) orient (PrimState m) m
zhpmv = hpmvAbstraction "zhpmv" cblas_zhpmv_safe cblas_zhpmv_unsafe withRStorable_

chpr :: PrimMonad m => HprFun Float (Complex Float) orient (PrimState m) m
chpr = hprAbstraction "chpr" cblas_chpr_safe cblas_chpr_unsafe (\x f -> f x)

zhpr :: PrimMonad m => HprFun Double (Complex Double) orient (PrimState m) m
zhpr = hprAbstraction "zhpr" cblas_zhpr_safe cblas_zhpr_unsafe (\x f -> f x)

chpr2 :: PrimMonad m => Hpr2Fun (Complex Float) orient (PrimState m) m
chpr2 = hpr2Abstraction "chpr2" cblas_chpr2_safe cblas_chpr2_unsafe withRStorable_

zhpr2 :: PrimMonad m => Hpr2Fun (Complex Double) orient (PrimState m) m
zhpr2 = hpr2Abstraction "zhpr2" cblas_zhpr2_safe cblas_zhpr2_unsafe withRStorable_

ssbmv :: PrimMonad m => SbmvFun Float orient (PrimState m) m
ssbmv = sbmvAbstraction "ssbmv" cblas_ssbmv_safe cblas_ssbmv_unsafe (\x f -> f x)

dsbmv :: PrimMonad m => SbmvFun Double orient (PrimState m) m
dsbmv = sbmvAbstraction "dsbmv" cblas_dsbmv_safe cblas_dsbmv_unsafe (\x f -> f x)

sspmv :: PrimMonad m => SpmvFun Float orient (PrimState m) m
sspmv = spmvAbstraction "sspmv" cblas_sspmv_safe cblas_sspmv_unsafe (\x f -> f x)

dspmv :: PrimMonad m => SpmvFun Double orient (PrimState m) m
dspmv = spmvAbstraction "dspmv" cblas_dspmv_safe cblas_dspmv_unsafe (\x f -> f x)

sspr :: PrimMonad m => SprFun Float orient (PrimState m) m
sspr = sprAbstraction "sspr" cblas_sspr_safe cblas_sspr_unsafe (\x f -> f x)

dspr :: PrimMonad m => SprFun Double orient (PrimState m) m
dspr = sprAbstraction "dspr" cblas_dspr_safe cblas_dspr_unsafe (\x f -> f x)

sspr2 :: PrimMonad m => Spr2Fun Float orient (PrimState m) m
sspr2 = spr2Abstraction "sspr2" cblas_sspr2_safe cblas_sspr2_unsafe (\x f -> f x)

dspr2 :: PrimMonad m => Spr2Fun Double orient (PrimState m) m
dspr2 = spr2Abstraction "dspr2" cblas_dspr2_safe cblas_dspr2_unsafe (\x f -> f x)

ssymv :: PrimMonad m => SymvFun Float orient (PrimState m) m
ssymv = symvAbstraction "ssymv" cblas_ssymv_safe cblas_ssymv_unsafe (\x f -> f x)

dsymv :: PrimMonad m => SymvFun Double orient (PrimState m) m
dsymv = symvAbstraction "dsymv" cblas_dsymv_safe cblas_dsymv_unsafe (\x f -> f x)

ssyr :: PrimMonad m => SyrFun Float orient (PrimState m) m
ssyr = syrAbstraction "ssyr" cblas_ssyr_safe cblas_ssyr_unsafe (\x f -> f x)

dsyr :: PrimMonad m => SyrFun Double orient (PrimState m) m
dsyr = syrAbstraction "dsyr" cblas_dsyr_safe cblas_dsyr_unsafe (\x f -> f x)

ssyr2 :: PrimMonad m => Syr2Fun Float orient (PrimState m) m
ssyr2 = syr2Abstraction "ssyr2" cblas_ssyr2_safe cblas_ssyr2_unsafe (\x f -> f x)

dsyr2 :: PrimMonad m => Syr2Fun Double orient (PrimState m) m
dsyr2 = syr2Abstraction "dsyr2" cblas_dsyr2_safe cblas_dsyr2_unsafe (\x f -> f x)

stbmv :: PrimMonad m => TbmvFun Float orient (PrimState m) m
stbmv = tbmvAbstraction "stbmv" cblas_stbmv_safe cblas_stbmv_unsafe

dtbmv :: PrimMonad m => TbmvFun Double orient (PrimState m) m
dtbmv = tbmvAbstraction "dtbmv" cblas_dtbmv_safe cblas_dtbmv_unsafe

ctbmv :: PrimMonad m => TbmvFun (Complex Float) orient (PrimState m) m
ctbmv = tbmvAbstraction "ctbmv" cblas_ctbmv_safe cblas_ctbmv_unsafe

ztbmv :: PrimMonad m => TbmvFun (Complex Double) orient (PrimState m) m
ztbmv = tbmvAbstraction "ztbmv" cblas_ztbmv_safe cblas_ztbmv_unsafe

stbsv :: PrimMonad m => TbsvFun Float orient (PrimState m) m
stbsv = tbsvAbstraction "stbsv" cblas_stbsv_safe cblas_stbsv_unsafe

dtbsv :: PrimMonad m => TbsvFun Double orient (PrimState m) m
dtbsv = tbsvAbstraction "dtbsv" cblas_dtbsv_safe cblas_dtbsv_unsafe

ctbsv :: PrimMonad m => TbsvFun (Complex Float) orient (PrimState m) m
ctbsv = tbsvAbstraction "ctbsv" cblas_ctbsv_safe cblas_ctbsv_unsafe

ztbsv :: PrimMonad m => TbsvFun (Complex Double) orient (PrimState m) m
ztbsv = tbsvAbstraction "ztbsv" cblas_ztbsv_safe cblas_ztbsv_unsafe

stpmv :: PrimMonad m => TpmvFun Float orient (PrimState m) m
stpmv = tpmvAbstraction "stpmv" cblas_stpmv_safe cblas_stpmv_unsafe

dtpmv :: PrimMonad m => TpmvFun Double orient (PrimState m) m
dtpmv = tpmvAbstraction "dtpmv" cblas_dtpmv_safe cblas_dtpmv_unsafe

ctpmv :: PrimMonad m => TpmvFun (Complex Float) orient (PrimState m) m
ctpmv = tpmvAbstraction "ctpmv" cblas_ctpmv_safe cblas_ctpmv_unsafe

ztpmv :: PrimMonad m => TpmvFun (Complex Double) orient (PrimState m) m
ztpmv = tpmvAbstraction "ztpmv" cblas_ztpmv_safe cblas_ztpmv_unsafe

stpsv :: PrimMonad m => TpsvFun Float orient (PrimState m) m
stpsv = tpsvAbstraction "stpsv" cblas_stpsv_safe cblas_stpsv_unsafe

dtpsv :: PrimMonad m => TpsvFun Double orient (PrimState m) m
dtpsv = tpsvAbstraction "dtpsv" cblas_dtpsv_safe cblas_dtpsv_unsafe

ctpsv :: PrimMonad m => TpsvFun (Complex Float) orient (PrimState m) m
ctpsv = tpsvAbstraction "ctpsv" cblas_ctpsv_safe cblas_ctpsv_unsafe

ztpsv :: PrimMonad m => TpsvFun (Complex Double) orient (PrimState m) m
ztpsv = tpsvAbstraction "ztpsv" cblas_ztpsv_safe cblas_ztpsv_unsafe

strmv :: PrimMonad m => TrmvFun Float orient (PrimState m) m
strmv = trmvAbstraction "strmv" cblas_strmv_safe cblas_strmv_unsafe

dtrmv :: PrimMonad m => TrmvFun Double orient (PrimState m) m
dtrmv = trmvAbstraction "dtrmv" cblas_dtrmv_safe cblas_dtrmv_unsafe

ctrmv :: PrimMonad m => TrmvFun (Complex Float) orient (PrimState m) m
ctrmv = trmvAbstraction "ctrmv" cblas_ctrmv_safe cblas_ctrmv_unsafe

ztrmv :: PrimMonad m => TrmvFun (Complex Double) orient (PrimState m) m
ztrmv = trmvAbstraction "ztrmv" cblas_ztrmv_safe cblas_ztrmv_unsafe

strsv :: PrimMonad m => TrsvFun Float orient (PrimState m) m
strsv = trsvAbstraction "strsv" cblas_strsv_safe cblas_strsv_unsafe

dtrsv :: PrimMonad m => TrsvFun Double orient (PrimState m) m
dtrsv = trsvAbstraction "dtrsv" cblas_dtrsv_safe cblas_dtrsv_unsafe

ctrsv :: PrimMonad m => TrsvFun (Complex Float) orient (PrimState m) m
ctrsv = trsvAbstraction "ctrsv" cblas_ctrsv_safe cblas_ctrsv_unsafe

ztrsv :: PrimMonad m => TrsvFun (Complex Double) orient (PrimState m) m
ztrsv = trsvAbstraction "ztrsv" cblas_ztrsv_safe cblas_ztrsv_unsafe

max_stars_repo_name: schnecki/hblas

max_issues_repo_name: schnecki/hblas

max_forks_repo_name: schnecki/hblas


text:
module TextStatistics (
    TextStatistics,
    wordCount,
    lineCount,
    averageWordLength,
    mostCommonLetter
) where

import Prelude hiding (filter)
import qualified Statistics as S
import Control.Applicative ((<$>), (<*>))
import qualified Data.HashMap.Strict as HM
import Data.Char

-- Statistics on Char stream

type TextStatistics a = S.Statistics Char a

-- predicates

isWordChar :: Char -> Bool
isWordChar = not <$> isSpace

isNewLine :: Char -> Bool
isNewLine c = c == '\n'

-- word count

data WordCountState = WordCountState {
    charCount :: Int,
    isInWord :: Bool
}

wordCount :: TextStatistics Int
wordCount = charCount <$> S.fold f (WordCountState 0 False)
    where
        f w c = WordCountState (if not (isInWord w) && isWordChar c then (charCount w) + 1 else (charCount w)) (isWordChar c)

-- line count. Assuming empty file has 1 line of 0 length

lineCount :: TextStatistics Int
lineCount = S.count isNewLine 1

-- average word length

wordCharCount :: TextStatistics Int
wordCharCount = S.count isWordChar 0

averageWordLength :: TextStatistics (Maybe Float)
averageWordLength = maybeDivide <$> wordCharCount <*> wordCount
    where
        maybeDivide a b = if b /= 0 then Just (fromIntegral a / fromIntegral b) else Nothing

-- most common character (if any) with number of usages

mostCommonLetter :: TextStatistics (Maybe Char)
mostCommonLetter = fmap fst <$> ((HM.foldlWithKey' maxChar Nothing) <$> lettersMap)
    where
        lettersMap :: TextStatistics (HM.HashMap Char Int)
        lettersMap = S.filter isWordChar $ S.fold (\hm c -> if HM.member c hm then HM.insertWith (+) c 1 hm else HM.insert c 1 hm) HM.empty
        maxChar Nothing c n = Just (c, n)
        maxChar (Just (c0, n0)) c n = if n > n0 then Just (c, n) else Just (c0, n0)

max_stars_repo_name: maciasoft/wc

max_issues_repo_name: maciasoft/wc

max_forks_repo_name: maciasoft/wc


text:
-- General utilities for working with neural networks.
--
-- Original author: David Banas <capn.freako@gmail.com>
-- Original date:   January 22, 2018
--
-- Copyright (c) 2018 David Banas; all rights reserved World wide.

{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_GHC -Wno-unused-top-binds #-}

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeSynonymInstances #-}

{-|
Module      : Haskell_ML.Util
Description : Provides certain general purpose utilities in the Haskell_ML package.
Copyright   : (c) David Banas, 2018
License     : BSD-3
Maintainer  : capn.freako@gmail.com
Stability   : experimental
Portability : ?
-}
module Utils.MLUtils where

import Prelude hiding (zipWith, zip)

import GHC.Generics (Par1(..),(:*:)(..),(:.:)(..))

import           Data.Finite
import           Data.Foldable
import           Data.Function                       (on)
import           Data.Key
import           Data.List                    hiding (zipWith, zip)
import           Data.Ord                            (comparing)
import           GHC.TypeLits
import qualified Data.Vector.Sized            as VS
import qualified Numeric.LinearAlgebra.Data   as LAD
import           Numeric.LinearAlgebra.Static hiding (mean)
import           Text.Printf

import           ConCat.Deep
import           ConCat.Misc                         (sqr, result, cond)
import           ConCat.Orphans                      (fstF, sndF)
import Utils.Utils

-- | Split a list of samples into training/testing sets.
--
-- The Finite given should be the percentage of samples desired for
-- training.
splitTrnTst :: Finite 101 -> [a] -> ([a],[a])
splitTrnTst n xs = splitAt (length xs * (fromInteger . getFinite) n `div` 100) xs

-- | Convert vector of Doubles from sized vector to hmatrix format.
toR :: KnownNat n => V n Double -> R n
toR = vector . VS.toList

-- | Calculate the classification accuracy, given:
--
--   - a functor of results vectors, and
--   - a functor of reference vectors.
classificationAccuracy :: (Foldable f, Keyed f, Eq (Key f), Foldable g, Zip g, Ord a)
                       => g (f a)  -- ^ Functor of result vectors.
                       -> g (f a)  -- ^ Functor of reference vectors.
                       -> Double
classificationAccuracy = (result.result) boolScore (zipWith ((==) `on` maxIndex))
 where boolScore = mean . fmap (cond 1 0)

-- | Mean value of a collection
mean :: (Foldable f, Fractional a) => f a -> a
mean = uncurry (/) . foldr (\e (s,c) -> (e+s,c+1)) (0,0)

-- | Prompt the user for input and read it.
prompt :: Read a => String -> IO a
prompt s = do
  putStrLn s
  readLn

-- | Pretty printer for values of type `R` n.
printVector :: (KnownNat n) => R n -> String
printVector v = foldl' (\ s x -> s ++ printf "%+6.4f  " x) "[ " ((LAD.toList . extract) v) ++ " ]"

-- | Pretty printer for values of type (`R` `m`, `R` `n`).
printVecPair :: (KnownNat m, KnownNat n) => (R m, R n) -> String
printVecPair (u, v) = "( " ++ printVector u ++ ", " ++ printVector v ++ " )"

-- | Show the learning evolution for the various parts/layers in a
-- list of HasLayers.
showPart :: (HasLayers p, RealFrac s, Ord s, PrintfArg s)
         => (p s -> [[s]])  -- ^ The "getter", which grabs the correct part of the parameters.
         -> String          -- ^ A label, which identifies the target of the getter.
         -> [p s]           -- ^ List of parameters, one per epoch.
         -> [String]
showPart f s ps =
  let xss = zip [1::Int,2..] $ (transpose . map f) ps
   in concatMap ( \ (i, xs) ->
                  [ "Average variance in layer " ++ show i ++ " " ++ s ++ ":"
                  , asciiPlot $ map (mean . map sqr) xs
                  ]
                ) xss

-- | Plot a list of Fractionals to an ASCII terminal.
asciiPlot :: (RealFrac a, Ord a, PrintfArg a) => [a] -> String
asciiPlot ys = unlines $
  zipWith (++)
    [ "        "
    , printf " %6.1e " x_max
    , "        "
    , "        "
    , "        "
    , "        "
    , "        "
    , "        "
    , "        "
    , "        "
    , "        "
    , printf " %6.1e " x_min
    , "        "
    ] $
    (:) "^" $ transpose (
    (:) "|||||||||||" $
    for xs $ \x ->
      valToStr $ (x - x_min) * 10 / x_range
    ) ++ ["|0" ++ concat [replicate 16 '_' ++ printf "%4d" (n * length ys `div` 3) | n <- [1..3]] ++ ">"]
      where valToStr  :: (RealFrac a, Ord a) => a -> String
            valToStr x = let i = round (10 - x)
                          in replicate i ' ' ++ "*" ++ replicate (10 - i) ' '
            x_min      = minimum xs
            x_max      = maximum xs
            x_range    = x_max - x_min
            xs         = takeEvery (length ys `div` 60) ys

-- | Return every Nth element of a list.
takeEvery :: Int -> [a] -> [a]
takeEvery n xs = case drop (n-1) xs of
                   (y:ys) -> y : takeEvery n ys
                   []     -> []

-- | Convenience function (= flip map).
for :: Functor f => f a -> (a -> b) -> f b
for = flip fmap

-- | Find the index of the maximum value in a keyed functor.
maxIndex :: (Foldable f, Keyed f, Ord a) => f a -> Key f
maxIndex = fst . maximumBy (comparing snd) . keyed

{--------------------------------------------------------------------
    Network structure
--------------------------------------------------------------------}

-- | A class of parameter types that can be split into a sequence of
-- layers of "weights" and "biases". The meaning of the two notions:
-- "weight" and "bias", are type specific. The only unifyng feature
-- of these, across types, is that they be fully resolved types.
class HasLayers p where
  getWeights :: p s -> [[s]]
  getBiases  :: p s -> [[s]]

instance (HasLayers f, HasLayers g) => HasLayers (g :*: f) where
  getWeights (g :*: f) = getWeights f ++ getWeights g
  getBiases  (g :*: f) = getBiases  f ++ getBiases  g

instance (Foldable a, Foldable b) => HasLayers (a --+ b) where
  getWeights (Comp1 gf) = [(concatMap (toList          . fstF) . toList) gf]
  getBiases  (Comp1 gf) = [(concatMap ((: []) . unPar1 . sndF) . toList) gf]

-- instance (HasLayers f, Foldable g) => HasLayers (g :.: f) where
--   getWeights (Comp1 g) = reverse $ foldl' (\ws -> (: ws) . concat . getWeights) [] g
--   getBiases  (Comp1 g) = reverse $ foldl' (\ws -> (: ws) . concat . getBiases)  [] g

-- instance Foldable f => HasLayers (f :*: Par1) where
--   getWeights (f :*: Par1 _) = [toList f]
--   getBiases  (_ :*: Par1 x) = [[x]]

-- instance Foldable f => HasLayers f where
--   getWeights f = [toList f]
--   getBiases  f = [[]]

max_stars_repo_name: faezs/opt-expect

max_issues_repo_name: faezs/opt-expect

max_forks_repo_name: faezs/opt-expect


text:
module Lib where

import Data.List ()
import Numeric.LinearAlgebra
  ( Linear (scale),
    Matrix,
    Vector,
    cmap,
    maxIndex,
    outer,
    (#>),
  )
import Text.CSV ()

{-
Sooooooooooooo
What do I need to do
  Need to make the output a vector
  make the weights a matrix
  make the learning stage be an outer product thats then
  subracted from the weights

  I need to implement sigmoid?
-}

predict ::
  Matrix Double -> -- curried - Weights
  Vector Double -> -- curried - Input vector
  Vector Double -- Guess
predict w x = cmap sigmoid w #> x

predictEpoch ::
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of input vectorors
  [Vector Double] -- A guess for each input in the epoch
predictEpoch w x = map (cmap sigmoid) (fmap (w #>) x)

sigmoid :: Double -> Double
sigmoid x = 1.0 / (1.0 + exp (negate x))

stocTrain ::
  Matrix Double -> -- w - Weights
  Vector Double -> -- x - Single input vectoror
  Vector Double -> -- l - label
  Matrix Double -- Trained Weights
stocTrain w x l = w - delta
  where
    error = predict w x - l
    sigDeriv = error * cmap (1 -) error
    delta = scale 0.001 $ sigDeriv `outer` x

trainEpoch ::
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of inputs
  [Vector Double] -> -- l - An epoch of labels
  Matrix Double -- Weights
trainEpoch w (x : xs) (l : ls) = trainEpoch (stocTrain w x l) xs ls
trainEpoch w [] [] = w

train ::
  Int -> -- t - number of epochs to train on
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of inputs
  [Vector Double] -> -- l - An epoch of labels
  Matrix Double -- Weights
train t w xs ls
  | t == 0 = w
  | otherwise = train (t - 1) (trainEpoch w xs ls) xs ls

successRate ::
  Matrix Double -> -- w - Weights
  [Vector Double] -> -- x - An epoch of input vectors
  [Vector Double] -> -- l - An epoch of labels
  Double -- A guess for each input in the epoch
successRate w x l = intDiv correct total
  where
    correct = sum (zipWith checkSame l (predictEpoch w x))
    total = length l

intDiv :: Int -> Int -> Double
intDiv n l = fromIntegral n / fromIntegral l

-- validate w x l = fmap maxIndex (predictEpoch w x)

checkSame :: Vector Double -> Vector Double -> Int
checkSame x y
  | maxIndex x == maxIndex y = 1
  | otherwise = 0

{-
so we need to:
stick the ending bias 1 onto each input array
turn the weights into a matrix
turn the input vectors in a list of vectors
turn the labels into a list of one hot encoded vectors
-}

-- csvToTrainingPred :: [[String]] -> [Double]
-- csvToTrainingPred s = predictEpoch trainedWeights scaledInputs
--   where
--     (labels, unscaledInputs) = labelInputSplit (map (map read) (init s))
--     scaledInputs = scaleInput unscaledInputs
--     trainedWeights = train 1 scaledInputs labels [0 ..]

-- scaleInput :: [[Double]] -> [[Double]]
-- scaleInput = map (map (/ 255))

-- labelInputSplit :: [[Double]] -> ([Double], [[Double]])
-- labelInputSplit xs = (head t, transpose (tail t))
--   where
--     t = transpose xs


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE Rank2Types #-}
{-# LANGUAGE FlexibleInstances    #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE IncoherentInstances  #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ForeignFunctionInterface #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}

module Data.Vector.Image.Color.YUV (
  YUV(..)
) where

import Prelude hiding(map)
import qualified Prelude as P
import qualified Data.List as L
import Data.Vector.Storable hiding(forM_,map,replicate,fromList,toList,(!?),(!))
import Data.Maybe
import Numeric.LinearAlgebra.Data hiding (fromList,toList,(!))
import Numeric.LinearAlgebra.HMatrix hiding (fromList,toList,(!))
import qualified Data.Vector.Storable as V
import Foreign
import GHC.Real
import Foreign.C.String
import Foreign.C.Types
import qualified Control.Monad as C
import Data.Vector.Image.Color.RGB

data YUV a = YUV {
   yuv_y :: a
 , yuv_u :: a
 , yuv_b :: a
} deriving (Show,Eq,Ord)

instance Functor YUV where
  fmap func (YUV r g b) = YUV (func r) (func g) (func b)

rgb2yuv :: (Num a,Integral a) => RGB a -> YUV a
rgb2yuv (RGB r g b) = YUV ((r+2*g+b)`div`4) (b-g) (r-g)

yuv2rgb :: (Num a,Integral a) => YUV a -> RGB a
yuv2rgb (YUV y u v) = RGB (v+g) g (u+g)
  where
    g = y - ((u + v)`div`4)

instance (Num a,Integral a) => ColorSpace YUV a where
  toRGB = yuv2rgb
  fromRGB = rgb2yuv
  toList (YUV r g b) = [r,g,b]
  fromList [r,g,b] = (YUV r g b)
  fromList val = error $ "List length must be 3 but " L.++  show (L.length val)
  dim _ = 3

max_stars_repo_name: junjihashimoto/img-vector

max_issues_repo_name: junjihashimoto/img-vector

max_forks_repo_name: junjihashimoto/img-vector


text:
{-# LANGUAGE FlexibleContexts #-}

module GT.Algebra.Matrix
  ( adjacencyMatrix,
    degreeMatrix,
    laplacianMatrix,
  )
where

import Data.Maybe (fromJust)
import Data.Sequence as Q (Seq, elemIndexL)
import Data.Tuple (swap)
import Data.Witherable (Witherable, mapMaybe)
import GT.Graph.Class
  ( Direction (Directed, Undirected),
    Graph (degree, direction, edgeMap, nodeCount, nodeIds, nodeMap),
    Unwrap (unwrap),
  )
import Numeric.LinearAlgebra.Data as LA
  ( Konst (konst),
    Matrix,
    R,
    accum,
    diagl,
  )

adjacencyMatrix :: Graph g n n' e e' d => g -> Matrix R
adjacencyMatrix g = LA.accum (konst 0 (n, n)) (+) as'
  where
    n = nodeCount g
    as' = case direction g of
      Directed -> as
      Undirected -> as ++ fmap (\(rc, a) -> (swap rc, a)) as
    is = nodeIds g
    as =
      edgeMap
        ( \e' ->
            let (i, j) = unwrap e'
                r = fromJust $ elemIndexL i is
                c = fromJust $ elemIndexL j is
             in ((r, c), 1)
        )
        g

degreeMatrix :: Graph g n n' e e' d => g -> Matrix R
degreeMatrix g = diagl $ nodeMap (\n' -> fromIntegral (fromJust $ degree (unwrap n') g)) g

laplacianMatrix :: Graph g n n' e e' d => g -> Matrix R
laplacianMatrix g = degreeMatrix g - adjacencyMatrix g

max_stars_repo_name: mskmysht/graphlib-hs

max_issues_repo_name: mskmysht/graphlib-hs

max_forks_repo_name: mskmysht/graphlib-hs