text:
#!/usr/bin/env stack
-- stack --resolver lts-15.04 runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
module Main where

import           Codec.Picture.Types
import           Control.Lens                  ()
import           Control.Monad
import           Data.Function
import           Data.List
import           Data.List.NonEmpty            (NonEmpty)
import qualified Data.List.NonEmpty            as NE
import           Data.Maybe
import           Data.Ratio
import qualified Data.Text                     as T
import           Data.Tuple
import qualified Data.Vector                   as V
import           Debug.Trace
import           Linear.Matrix                 hiding (trace)
import           Linear.Metric
import           Linear.V2
import           Linear.V3
import           Linear.Vector
import           Numeric.LinearAlgebra         hiding (polar, scale, (<>))
import qualified Numeric.LinearAlgebra         as Matrix
import           Numeric.LinearAlgebra.HMatrix hiding (polar, scale, (<>))
import           Reanimate
import           Reanimate.Animation
import           Reanimate.Math.Balloon
import           Reanimate.Math.Common
import           Reanimate.Math.Triangulate
import           Reanimate.Math.Polygon
import           Reanimate.Math.EarClip
import           Reanimate.Math.SSSP
import           Reanimate.Math.Render
import           Reanimate.Math.Visibility
import           Reanimate.Math.Compatible
import           Reanimate.Morph.Common
import           Reanimate.PolyShape           (svgToPolygons)

p1 = pScale 2 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
-- p1 = pSetOffset (addPoints 2 shape13) 0
-- p2 = pScale 0.5 $ centerPolygon shape20
-- p1 = centerPolygon shape2
p2 = pScale 2 $ pAtCenter $ pAddPoints 0 (pSetOffset shape14 0)

main :: IO ()
main = reanimate $ sceneAnimation $ do
  bg <- newSpriteSVG $ mkBackground "black"
  spriteZ bg (-1)
  -- play $ drawVisibleFrom triangle
  -- play $ drawVisibleFrom shape1
  -- play $ drawVisibleFrom shape2
  -- play $ drawVisibleFrom shape3
  -- play $ drawVisibleFrom shape4
  -- play $ drawVisibleFrom shape5
  -- play $ drawVisibleFrom shape6
  -- play $ drawOverlap shape20
  -- play $ drawSSSP triangle naive
  -- play $ drawSSSPFast shape5
  -- fork $ play $ mapA (translate (-3) 0) $ drawSSSPVisibilityFast shape2
  -- fork $ play $ mapA (translate (3) 0) $ drawSSSPVisibilityFast shape7
  -- play $ drawOverlap $ fst $ split1Link (fst (split1Link shape7 3 7)) 0 2
  -- fork $ play $ mapA (translate (-4) 0) $ drawSSSPVisibilityFast $ pSetOffset (addPoints 0 shape14) 0
  -- fork $ play $ mapA (translate (4) 0) $ drawSSSPVisibilityFast $ pSetOffset (addPoints 2 shape13) 0
  fork $ play $ staticFrame 1 $
    translate (4) 0 $ mkGroup
    [ withFillColor "grey" $ polygonShape p1
    , polygonNumDots p1 ]
  fork $ play $ staticFrame 1 $
    let -- pOrigin = (addPoints 10 $ pSetOffset shape14 0)
        pOrigin = pScale 2 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
        --pOrigin = (addPoints 2 shape13)
        p1 = pScale 1 $ pSetOffset pOrigin 0
        p2 = snd $ split1Link p1 1 11 0
        p3 = fst $ split1Link p2 0 2 1
        -- p4 = fst $ split2Link p3 2 5
        p4 = snd $ split1Link p3 0 9 1
        p5 = fst $ split1Link p4 1 3 0
        p6 = fst $ split2Link p5 0 2
        p7 = fst $ split1Link p6 1 3 0
        p8 = fst $ split2Link p7 0 2
        p9 = fst $ split1Link p7 1 3 0
        p = p4
    in mkGroup
    [ withFillColor "grey" $ polygonShape p
    , polygonNumDots p ]
  -- newSpriteSVG $
  --   let p1 = pSetOffset shape14 0
  --       V2 x y = realToFrac <$> steiner2Link p1 2 6
  --   in translate (-4) 0 $
  --     translate x y $ withFillColor "red" $
  --     mkCircle 0.1
  -- newSpriteSVG $
  --   let p1 = pSetOffset shape14 0
  --   in translate (-4) 0 $
  --     mkGroup
  --     [ mkGroup
  --       [ withStrokeColor "purple" $
  --         mkLinePath [(x1,y1),(x2,y2)]
  --       , translate x2 y2 $ withFillColor "red" $ mkCircle 0.1 ]
  --     | (eA, eB) <- take 1 $ steiner2Edges p1 4 11
  --     , let V2 x1 y1 = realToFrac <$> eA
  --           V2 x2 y2 = realToFrac <$> eB
  --     ]
  -- play $ staticFrame 1 $
  --   let p1 = pSetOffset shape14 0
  --   in mkGroup
  --   [ withFillColor "grey" $ polygonShape p1
  --   , polygonNumDots p1
  --   , drawWindowOverlap p1 1 6
  --   -- , withStrokeColor "red" $ drawWindow (pSetOffset p1 2)
  --   -- , withStrokeColor "blue" $ drawWindow (pSetOffset p1 6)
  --   ]
  -- play $ drawCompatible (pSetOffset (addPoints 2 shape13) 0) (pSetOffset shape14 0)

  -- play $ drawSSSP shape2 naive
  -- play $ drawSSSP shape3 naive
  -- play $ drawSSSP shape4 naive
  -- play $ drawSSSP shape5 naive
  -- play $ drawSSSP (pScale 0.5 $ winding 10) (\p -> sssp p (dual (earClip p)))
  -- play $ drawSSSP shape1 (\p -> sssp p (dual (earClip p)))
  -- play $ drawTriangulation (pCycle shape5 0.2910962834555265) earClip'
  -- play $ drawTriangulation shape5 earClip'
  -- play $ mkAnimation 1 $ \t ->
  --   let p = shape4
  --   in polygonNumDots (pCycle p t)
  -- play $ setDuration 20 $ drawSSSPVisibility $ pScale 1 $ shape7
  -- let shapeI = head $ svgToPolygons 0.1 $ scale 8 $ center $ latex "I"
  -- play $ animate $ \_ -> withFillColor "grey" $ polygonNumDots shapeI
  -- play $ drawTriangulation (pScale 0.5 $ winding 10) earClip'
  -- play $ staticFrame 1 $
  --   mkGroup
  --   [ withFillColor "grey" $ polygonShape shape13
  --   , withFillColor "grey" $ polygonDots shape13 ]
  -- let p = balloonP origin
  --     origin = centerPolygon $ pScale 1 $ shiftLongestDiameter shapeI
  --     mkB = balloon (scale 8 $ center $ latex "C")
  --     inf = unsafeSVGToPolygon 0.1 $ (scale 8 $ center $ latex "$\\infty$")
  --     cShape = shiftLongestDiameter $ unsafeSVGToPolygon 0.01 $ (scale 8 $ center $ latex "C")
  -- _ <- newSpriteSVG $
  --   translate (0) (3) $ withFillColor "white" $
  --   center $ latex $ T.pack $ show (isSimple inf)
  -- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t ->
  --   let inflated = p t in
  --   translate (0) (0) $ withFillColor "white" $ withStrokeColor "red" $
  --   withStrokeWidth (defaultStrokeWidth*0) $
    -- polygonShape inf
    -- renderTriangulation inf (earClip inf)
    -- mkB t
    -- scale 3 $ -- translate (-2.2) (-1) $
    -- mkGroup
    -- [ mkGroup []
    -- , withStrokeWidth (defaultStrokeWidth*0.05) $
    --   renderDual (pRing shape22) (dual (polygonOffset shape22) $ polygonTriangulation shape22)
    -- , polygonNumDots shape22
    --   -- renderTriangulation cShape (polygonTriangulation cShape)
    -- ]
    -- mkGroup
    -- [ -- translate (-2) 0 $ withFillColor "white" $ polygonShape $ balloonP (min 1 $ t+0.1) origin
    --   if False then mkGroup [] else translate (0) 0 $ mkGroup
    --   [ withFillColor "white" $ polygonShape inflated
    --   -- , polygonNumDots inflated
    --   ]
    -- -- , translate (2) 0 $ mkGroup
    -- --   [ withFillColor "grey" $ polygonShape origin
    -- --   , polygonNumDots origin
    -- --   ]
    -- ]
  -- play $ mapA (withStrokeWidth (defaultStrokeWidth*0.2)) $ drawTriangulation cShape earClip'
  -- play $ staticFrame 1 $ renderTriangulation shape3 earClip
  -- play $ staticFrame 1 $ renderTriangulation shape4 earClip
  -- play $ staticFrame 1 $ renderTriangulation shape5 earClip
  -- play $ staticFrame 1 $ renderTriangulation shape6 earClip
  -- let p = deoverlapPolygon shape23
  -- newSpriteSVG $ translate (-3) 0 $ scale 2 $ center $ mkGroup
  --   [ withFillColor "grey" $ polygonShape p
  --   , polygonNumDots p
  --   ]
  -- wait 1
  -- play $ mapA (scale 2 . withStrokeWidth (defaultStrokeWidth*0.2)) $
  --   drawTriangulation p (earClip')
  return ()

drawVisibility :: Polygon -> Animation
drawVisibility p' = mkAnimation 5 $ \t ->
  let p = pCycle p' (t::Double) in
  centerUsing (polygonShape p) $
  mkGroup
  [ withFillColor "grey" $ polygonShape p
  , withFillColor "grey" $ polygonDots p
  , withFillColor "white" $ mkLinePathClosed
    [ (x,y) | V2 x y <- visibility (map (fmap realToFrac) $ V.toList $ polygonPoints p) ]
  , let V2 x y = fmap realToFrac $ pAccess p 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.1
  -- , withFillColor "blue" $ latex $ T.pack $ show (t)
  ]

drawSSSPVisibility :: Polygon -> Animation
drawSSSPVisibility p' = mkAnimation 5 $ \t ->
  let p = pSetOffset p' (round $ t*(fromIntegral $ pSize p'-1))
      vis = ssspVisibility p in
  centerUsing (polygonShape p) $
  mkGroup
  [ withFillColor "grey" $ polygonShape p
  -- , withFillColor "grey" $ polygonDots p
  -- , withFillColor "white" $ polygonShape vis
  , let V2 x y = fmap realToFrac $ pAccess p 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.1
  -- , withFillColor "blue" $ latex $ T.pack $ show (t)
  ]

drawSSSPVisibilityFast :: Polygon -> Animation
drawSSSPVisibilityFast p' = mkAnimation 5 $ \t ->
  let root = min (pSize p-1) $ (floor $ t*(fromIntegral $ pSize p))
      p = pSetOffset p' root
      vis = ssspVisibility p in
  -- centerUsing (polygonShape p) $
  mkGroup
  [ withFillColor "grey" $ polygonShape p
  , withFillColor "white" $ polygonShape vis
  , withFillColor "grey" $ polygonNumDots p
  , let V2 x y = fmap realToFrac $ pAccess p 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.09
  -- , withFillColor "blue" $ latex $ T.pack $ show (t)
  ]

drawCompatible :: Polygon -> Polygon -> Animation
drawCompatible a b = sceneAnimation $ do
  newSpriteSVG $ translate (-3) 0 $ mkGroup
    [ withFillColor "grey" $ polygonShape a
    , withFillColor "grey" $ polygonNumDots a
    ]
  newSpriteSVG $ translate (3) 0 $ mkGroup
    [ withFillColor "grey" $ polygonShape b
    , withFillColor "grey" $ polygonNumDots b
    ]
  let compat = compatiblyTriangulateP a b
  forM_ compat $ \(l, r) -> do
    fork $ play $ staticFrame 1 $
      translate (-3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
      withFillOpacity 0 $ polygonShape l
    fork $ play $ staticFrame 1 $
      translate (3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
      withFillOpacity 0 $ polygonShape r
  -- forM_ compat $ \(l, r) -> waitOn $ do
  --   fork $ play $ staticFrame 1 $
  --     translate (-3) 0 $
  --     withFillColor "white" $ polygonShape l
  --   fork $ play $ staticFrame 1 $
  --     translate (3) 0 $
  --     withFillColor "white" $ polygonShape r

drawOverlap :: Polygon -> Animation
drawOverlap p' = mkAnimation 5 $ \t ->
  let p = pCycle p' (t::Double)
      vis = ssspVisibility p
      vis' = ssspVisibility p'
      mWins = ssspWindows p
      oWins = ssspWindows p'
      -- (left, right) = split1Link p' 0 3
      -- (left, right) = split2Link p' 0 2
      -- (left, right) = splitNLink p' 0 [(TwoLink,2),(OneLink,3)]
      -- (left, right) = splitNLink p' 0 [(OneLink,5),(TwoLink,3)]
      -- sPoly = if t < 0.5 then left else right
      in
  centerUsing (polygonShape p) $
  mkGroup
  [ withFillColor "grey" $ polygonShape p
  -- , withFillColor "grey" $ polygonDots p
  , withFillOpacity 0.5 $ withFillColor "lightgreen" $ polygonShape vis
  -- , withFillOpacity 0.5 $ withFillColor "cyan" $ polygonShape vis'
  , let V2 x y = fmap realToFrac $ pAccess p 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.1
  , let V2 x y = fmap realToFrac $ pAccess p' 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.1
  -- , withFillColor "blue" $ latex $ T.pack $ show (t)
  , mkGroup
    [ withStrokeColor "red" $
      mkLine (x1,y1) (x2,y2)
    | (a,b) <- mWins
    , let V2 x1 y1 = realToFrac <$> a
          V2 x2 y2 = realToFrac <$> b
    ]
  -- , mkGroup
  --   [ withStrokeColor "blue" $
  --     mkLine (x1,y1) (x2,y2)
  --   | (a,b) <- oWins
  --   , let V2 x1 y1 = realToFrac <$> a
  --         V2 x2 y2 = realToFrac <$> b
  --   ]
  , mkGroup
    [ withStrokeColor "green" $mkGroup
      [ mkLine (x1,y1) (x2,y2)
      , mkLine (x1',y1') (x2',y2') ]
    | (a,b) <- mWins
    , (i,j) <- oWins
    , a == i || a == j || b == i || b == j
    , not (sort [a,b] == sort [i,j])
    , let V2 x1 y1 = realToFrac <$> a
          V2 x2 y2 = realToFrac <$> b
          V2 x1' y1' = realToFrac <$> i
          V2 x2' y2' = realToFrac <$> j
    ]
  -- , mkGroup
  --   [ let V2 x y = realToFrac <$> link in
  --     translate x y $
  --     withFillColor "red" $
  --     mkCircle 0.1
  --   | link <- [steiner2Link p' 0 2, steiner2Link p' 5 3] ]
  -- , withFillColor "blue" $ polygonShape sPoly
  ]

drawWindow :: Polygon -> SVG
drawWindow p =
  let mWins = ssspWindows p
      in
  mkGroup
  [ mkGroup
    [ mkLine (x1,y1) (x2,y2)
    | (a,b) <- mWins
    , let V2 x1 y1 = realToFrac <$> a
          V2 x2 y2 = realToFrac <$> b
    ]
  ]

drawWindowOverlap :: Polygon -> Int -> Int -> SVG
drawWindowOverlap p a b =
  let aWins = ssspWindows (pAdjustOffset p a)
      bWins = ssspWindows (pAdjustOffset p b)
      in
  mkGroup
  [ mkGroup $
    [ withStrokeColor "green" $mkGroup
      [ mkLine (x1,y1) (x2,y2)
      , mkLine (x1',y1') (x2',y2') ]
    | (a,b) <- aWins
    , (i,j) <- bWins
    , a == i || a == j || b == i || b == j
    , not (sort [a,b] == sort [i,j])
    , let V2 x1 y1 = realToFrac <$> a
          V2 x2 y2 = realToFrac <$> b
          V2 x1' y1' = realToFrac <$> i
          V2 x2' y2' = realToFrac <$> j
    ]
  ]

drawSSSP :: Polygon -> (Polygon -> SSSP) -> Animation
drawSSSP p gen = mkAnimation 5 $ \t -> centerUsing outline $
  let p' = pCycles p !! (round $ t*(fromIntegral $ pSize p-1)) in
  mkGroup
  [ outline
  , renderSSSP p' (gen p')
  -- , let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
  --   translate x y $ withFillColor "red" $ mkCircle 0.1
  , withFillColor "grey" $ polygonNumDots $ p'
  ]
  where
    outline =
      withFillColor "grey" $ mkLinePathClosed
        [ (x,y) | V2 x y <- map (fmap realToFrac) (V.toList (polygonPoints p)  ++ [pAccess p 0]) ]

{-# INLINE drawSSSPFast #-}
drawSSSPFast :: Polygon -> Animation
drawSSSPFast p = mkAnimation 5 $ \t -> centerUsing outline $
  let root = (round $ t*(fromIntegral $ pSize p-1))
      d = dual root triangulation
      sTree = sssp (pRing p) d in
  mkGroup
  [ outline
  , renderSSSP p sTree
  -- , let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
  --   translate x y $ withFillColor "red" $ mkCircle 0.1
  , withFillColor "grey" $ polygonNumDots $ p
  ]
  where
    triangulation = earClip $ pRing p
    outline =
      withFillColor "grey" $ mkLinePathClosed
        [ (x,y) | V2 x y <- map (fmap realToFrac) (V.toList (polygonPoints p) ++ [pAccess p 0]) ]

drawVisibleFrom :: Polygon -> Animation
drawVisibleFrom p = mkAnimation 5 $ \t -> centerUsing (polygonShape p) $ mkGroup
  [ withFillColor "grey" $ polygonShape p
  , withFillColor "grey" $ polygonDots p
  , renderVisibleFrom (pCycle p t)
  , let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
    translate x y $ withFillColor "red" $ mkCircle 0.1
  ]

renderVisibleFrom :: Polygon -> SVG
renderVisibleFrom p = withStrokeColor "white" $ withFillColor "white" $ mkGroup
  [ mkGroup
    [ mkLine (ax,ay) (bx,by)
    , translate bx by $ mkCircle 0.1 ]
  | i <- visibilityArray (pRing p) V.! 0
  , let V2 ax ay = fmap realToFrac $ pAccess p 0
        V2 bx by = fmap realToFrac $ pAccess p i ]

drawTriangulation :: Polygon -> (Ring Rational -> [Triangulation]) -> Animation
drawTriangulation p gen = sceneAnimation $ do
  forM_ (gen $ pRing p) $ \t -> play $ staticFrame 1 $ renderTriangulation p t

renderDual :: Ring Rational -> Dual -> SVG
renderDual ring d = case d of
    Dual (a,b,c) l r -> mkGroup
      [ withFillColor "blue" $ mkTrig a b c
      , worker c a l
      , worker b c r
      ]
  where
    mkTrig a b c =
      let V2 x1 y1 = realToFrac <$> ringAccess ring a
          V2 x2 y2 = realToFrac <$> ringAccess ring b
          V2 x3 y3 = realToFrac <$> ringAccess ring c
      in mkLinePathClosed [ (x1, y1), (x2,y2), (x3,y3) ]
    worker p1 p2 EmptyDual = mkGroup []
    worker p1 p2 (NodeDual x l r) = mkGroup
      [ mkTrig p1 p2 x
      , worker x p2 l
      , worker p1 x r
      ]
-- data Dual = Dual (Int,Int,Int) -- (a,b,c)
--                   DualTree -- borders ca
--                   DualTree -- borders bc
--   deriving (Show)
--
-- data DualTree
--   = EmptyDual
--   | NodeDual Int -- axb triangle, a and b are from parent.
--       DualTree -- borders xb
--       DualTree -- borders ax
--   deriving (Show)

max_stars_repo_name: TristanCacqueray/reanimate

max_issues_repo_name: TristanCacqueray/reanimate

max_forks_repo_name: TristanCacqueray/reanimate


text:
module Math.TestDimensionsPCA where

import Universum hiding (All, Vector, transpose)

import GPLVM.Types

import Control.Lens (makeLenses)
import Data.Array.Repa
import Data.Array.Repa.Repr.ForeignPtr
import Numeric.LinearAlgebra.Repa hiding (Matrix, Vector)

import Math.Matrix

data PCA = PCA
  { inputData    :: Matrix D Double
  , covariance   :: Covariance Double
  , eigenVectors :: EigenVectors Double
  , eigenValues  :: Vector F Double
  , finalData    :: Matrix D Double
  , restoredData :: Matrix D Double
  }

makePCA
    :: Int
    -> Matrix D Double
    -> PCA
makePCA desiredDimensions input =
  let inputData = input
      (Z :. yInp :. _) = extent input                                          -- get dimension
      (meanVec, covariance) = meanCovS input
      meanMatrix = DimMatrix $ toMatrix meanVec yInp
      adjustInput = DimMatrix input -^^ meanMatrix
      (eigenValues', eigenVec) = eigSHM covariance
      eigenValue = runVector eigenValues'
      eigenVecD@(DimMatrix (ADelayed (Z :. y :. x) f)) = transposeM eigenVec              -- transposing is almost free
      eigenVectors' =                                                        -- leave only n desired eigenvectors
        if desiredDimensions >= y
        then eigenVecD
        else DimMatrix $ ADelayed (Z :. desiredDimensions :. x) f
      eigenVectors = getInternal $ transposeM eigenVectors'                       -- colunmns are eigenvectors
      finalData' = eigenVectors' `mulPM` transposeM adjustInput -- data in new eigenvectors space
      finalData = getInternal . transposeM $ finalData'
      restoredDataWOMean = transposeM $ pinvSM eigenVectors' `mulPM` finalData'  -- restore to the original space
      restoredData = getInternal $ restoredDataWOMean +^^ meanMatrix              -- add mean
  in PCA{..}

ext: hs

max_stars_repo_name: serokell/GPLVMHaskell

max_issues_repo_name: serokell/GPLVMHaskell

max_forks_repo_name: serokell/GPLVMHaskell


text:
module HLearn.Models.Distributions.Multivariate.MultiNormalFast
    ( MultiNormal
    , covar
    )
    where

import Control.DeepSeq
import qualified Data.Foldable as F
import qualified Data.Semigroup as SG

import Foreign.Storable
import Numeric.LinearAlgebra hiding ((<>))
import qualified Numeric.LinearAlgebra as LA

import HLearn.Algebra
import HLearn.Models.Distributions.Common

-------------------------------------------------------------------------------
-- data types

type MultiNormal dp = AddUnit1 MultiNormal' dp

data MultiNormal' dp = MultiNormal'
    { raw       :: !(RawMoments dp)
    , meanV     :: Vector (Ring dp)
    , covarM    :: Matrix (Ring dp)
    , invcovarM :: Matrix (Ring dp)
    , pdfdenom  :: Ring dp
    }

deriving instance (Element (Ring dp), Read (Ring dp)) => Read (MultiNormal' dp)
deriving instance (Element (Ring dp), Show (Ring dp)) => Show (MultiNormal' dp)

instance (NFData (Ring dp), Storable (Ring dp)) => NFData (MultiNormal' dp) where
    rnf m = deepseq (raw m)
          $ deepseq (meanV m)
          $ deepseq (covarM m)
          $ deepseq (invcovarM m)
          $ rnf (pdfdenom m)

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

data RawMoments dp = RawMoments
    { m0 :: !(Ring dp)
    , m1 :: !(Vector (Ring dp))
    , m2 :: !(Matrix (Ring dp))
    }

deriving instance (Element (Ring dp), Read (Ring dp)) => Read (RawMoments dp)
deriving instance (Element (Ring dp), Show (Ring dp)) => Show (RawMoments dp)

instance NFData (RawMoments dp) where
    rnf rm = seq (m0 rm)
           $ seq (m1 rm)
           $ seq (m2 rm) ()

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

covar :: MultiNormal dp -> Matrix (Ring dp)
covar (UnitLift1 m) = covarM m

mkMultiNormal' ::
    ( Storable (Ring dp)
    , Element (Ring dp)
    , Field (Ring dp)
    , Floating (Ring dp)
    ) => RawMoments dp -> MultiNormal' dp
mkMultiNormal' raw = MultiNormal'
    { raw = raw
    , meanV = meanV'
    , covarM = covarM'
    , invcovarM = inv covarM'
    , pdfdenom = 1/(sqrt $ det covarM' * (2*pi)^d)
    }
    where
        n = m0 raw
        d = dim $ m1 raw
        meanV' = mapVector (/n) (m1 raw)
        covarM' = mapMatrixWithIndex (\(i,j) mij -> (mij/n-(meanV' @> j)*(meanV' @> i))) (m2 raw)

-------------------------------------------------------------------------------
-- algebra

instance
    ( Container Matrix (Ring dp)
    , Container Vector (Ring dp)
    , Field (Ring dp)
    , Floating (Ring dp)
    ) => SG.Semigroup (MultiNormal' dp)
        where
    d1 <> d2 = mkMultiNormal' $ RawMoments
        { m0 = m0 (raw d1) + m0 (raw d2)
        , m1 = m1 (raw d1) `LA.add` m1 (raw d2)
        , m2 = m2 (raw d1) `LA.add` m2 (raw d2)
        }

-------------------------------------------------------------------------------
-- training

instance Num a => HasRing [a] where
    type Ring [a] = a

instance
    ( Element r
    , Container Vector r
    , Field r
    , Floating r
    , Ring (dp r) ~ r
    , F.Foldable dp
    ) => HomTrainer (MultiNormal (dp r))
        where
    type Datapoint (MultiNormal (dp r)) = dp r

    train1dp dp = UnitLift1 . mkMultiNormal' $ RawMoments
        { m0 = 1
        , m1 = fromList $ F.toList dp
        , m2 = n><n $ [i*j | i <- F.toList dp, j <- F.toList dp]
        }
        where
            n = length $ F.toList dp

-------------------------------------------------------------------------------
-- distribution

instance Probabilistic (MultiNormal dp) where
    type Probability (MultiNormal dp) = Ring dp

instance
    ( Element r
    , Container Vector r
    , Floating r
    , LA.Product r
    , Ring (dp r) ~ r
    , F.Foldable dp
    ) => PDF (MultiNormal (dp r))
        where
    pdf (UnitLift1 dist) dp = pdfdenom dist * (exp $ (-1/2) * top)
        where
            x = fromList $ F.toList dp
            top=minElement $ ((asRow  x) `sub` (asRow $ meanV dist))
                       LA.<> (invcovarM dist)
                       LA.<> ((asColumn x) `sub` (asColumn $ meanV dist))

max_stars_repo_name: silky/HLearn

max_issues_repo_name: silky/HLearn

max_forks_repo_name: silky/HLearn


text:
{-# LANGUAGE RecordWildCards #-}

module Statistics.Block
    ( BlockHeader (..)
    , blockHeadersF
    , blockChain
    , blockChainF
    , txBlocksF
    , inBlockChainF
    , txCntInChainF
    , TxFate (..)
    , txFateF
    ) where

import           Control.Foldl (Fold (..), fold)
import qualified Data.Map.Lazy as ML
import qualified Data.Map.Strict as MS
import           Data.Maybe (fromJust, isJust)
import qualified Data.Set as S
import qualified Data.Text as T
import           Data.Time.Units (Microsecond)

import           JSONLog (IndexedJLTimedEvent (..))
import           Pos.Infra.Util.JsonLog.Events (JLBlock (..), JLEvent (..))
import           Prelude (id)
import           Statistics.Tx (txFirstReceivedF)
import           Types
import           Universum hiding (fold)

data BlockHeader = BlockHeader
    { bhNode      :: !NodeId
    , bhTimestamp :: !Timestamp
    , bhHash      :: !BlockHash
    , bhPrevBlock :: !BlockHash
    , bhSlot      :: !Slot
    , bhTxCnt     :: !Int
    } deriving Show

type SMap k a = MS.Map k a
type LMap k a = ML.Map k a

blockHeadersF :: Fold IndexedJLTimedEvent (SMap BlockHash BlockHeader)
blockHeadersF = Fold step MS.empty id
  where
    step :: SMap Text BlockHeader -> IndexedJLTimedEvent -> SMap Text BlockHeader
    step m IndexedJLTimedEvent{..} = case ijlEvent of
        JLCreatedBlock JLBlock{..} ->
            let bh = BlockHeader
                    { bhNode      = ijlNode
                    , bhTimestamp = ijlTimestamp
                    , bhHash      = jlHash
                    , bhPrevBlock = jlPrevBlock
                    , bhSlot      = jlSlot
                    , bhTxCnt     = length jlTxs
                    }
            in MS.insert jlHash bh m

        _                          -> m

blockChain :: SMap BlockHash BlockHeader -> Set BlockHash
blockChain m = S.fromList $ maybe [] id $ (fmap fst . uncons) $ sortByLengthDesc [getLongestChain h | h <- allHashes]
  where
    allHashes :: [BlockHash]
    allHashes = MS.keys m

    sortByLengthDesc :: [[a]] -> [[a]]
    sortByLengthDesc = sortBy (compare `on` (negate . length))

    successors :: SMap BlockHash [BlockHash]
    successors = foldl' f MS.empty allHashes

    getSuccessors :: BlockHash -> [BlockHash]
    getSuccessors h = MS.findWithDefault [] h successors

    f :: SMap BlockHash [BlockHash] -> BlockHash -> SMap BlockHash [BlockHash]
    f s h = let BlockHeader{..} = m MS.! h
            in  MS.alter (Just . maybe [bhHash] (bhHash :)) bhPrevBlock s

    longestChains :: LMap Text [Text]
    longestChains = ML.fromList [(h, getLongestChain h) | h <- allHashes]

    getLongestChain :: Text -> [Text]
    getLongestChain h = case sortByLengthDesc $ map (longestChains ML.!) $ getSuccessors h of
        []      -> [h]
        (c : _) -> h : c

blockChainF :: Fold IndexedJLTimedEvent (Set BlockHash)
blockChainF = blockChain <$> blockHeadersF

txBlocksF :: Fold IndexedJLTimedEvent (SMap TxHash [(Timestamp, BlockHash)])
txBlocksF = Fold step MS.empty id
  where
    step :: SMap Text [(Microsecond, Text)]
         -> IndexedJLTimedEvent
         -> SMap Text [(Microsecond, Text)]
    step m IndexedJLTimedEvent{..} = case ijlEvent of
        JLCreatedBlock JLBlock{..} -> foldl' (f ijlTimestamp jlHash) m [T.take 16 x | x <- jlTxs]
        _                          -> m

    f :: Timestamp
      -> Text
      -> SMap Text [(Timestamp, Text)]
      -> Text
      -> SMap Text [(Timestamp, Text)]
    f ts h m tx = let y = (ts, h)
                  in  MS.alter (Just . maybe [y] (y :)) tx m

inBlockChainF :: Fold IndexedJLTimedEvent (SMap TxHash Timestamp)
inBlockChainF = f <$> txFateF
  where
    f :: SMap TxHash TxFate -> SMap TxHash Timestamp
    f = MS.map fromJust . MS.filter isJust .  MS.map txInBlockChain

txCntInChainF :: Fold IndexedJLTimedEvent [(NodeId, Timestamp, Int)]
txCntInChainF = f <$> blockHeadersF <*> blockChainF
  where
    f :: Map BlockHash BlockHeader -> Set BlockHash -> [(NodeId, Timestamp, Int)]
    f m cs = [(bhNode, bhTimestamp, bhTxCnt) | (_, BlockHeader{..}) <- MS.toList m, bhHash `S.member` cs]

data TxFate =
      InNoBlock
    | InBlockChain !Timestamp !BlockHash !(Set (Timestamp, BlockHash))
    | InFork !(Set (Timestamp, BlockHash))
    deriving Show

txInBlockChain :: TxFate -> Maybe Timestamp
txInBlockChain InNoBlock             = Nothing
txInBlockChain (InBlockChain ts _ _) = Just ts
txInBlockChain (InFork _)            = Nothing

txFateF :: Fold IndexedJLTimedEvent (SMap TxHash TxFate)
txFateF = f <$> txFirstReceivedF <*> txBlocksF <*> blockChainF
  where
    f :: SMap TxHash Timestamp -> SMap TxHash [(Timestamp, BlockHash)] -> Set BlockHash -> SMap TxHash TxFate
    f received blocks chain = MS.fromList [(tx, fate tx) | tx <- MS.keys received]
      where
        fate :: TxHash -> TxFate
        fate tx = case MS.lookup tx blocks of
            Nothing -> InNoBlock
            Just xs -> fold (Fold step (Nothing, S.empty) extract) xs

        step :: (Maybe (Timestamp, BlockHash), Set (Timestamp, BlockHash))
             -> (Timestamp, BlockHash)
             -> (Maybe (Timestamp, BlockHash), Set (Timestamp, BlockHash))
        step (Nothing, s) (ts, h)
            | S.member h chain = (Just (ts, h), s)
            | otherwise        = (Nothing     , S.insert (ts, h) s)
        step (p      , s) q    = (p           , S.insert q       s)

        extract :: (Maybe (Timestamp, BlockHash), Set (Timestamp, BlockHash)) -> TxFate
        extract (Nothing     , s) = InFork s
        extract (Just (ts, h), s) = InBlockChain ts h s

max_stars_repo_name: Davdegreat/cardano-sl

max_issues_repo_name: Davdegreat/cardano-sl

max_forks_repo_name: Davdegreat/cardano-sl


text:
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE InstanceSigs #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Commons where

import Control.Applicative (liftA2)
import Control.Monad (foldM)
import Data.Array hiding (range)
import Data.Complex (Complex (..), imagPart, realPart)
import Data.Function ((&))
import Data.List (intercalate, sort)
import qualified Data.Map.Strict as Map
import Data.Set (fromList, toList)
import HashedExpression.Internal
import HashedExpression.Internal.Base
import HashedExpression.Internal.Builder
import HashedExpression.Internal.MonadExpression
import HashedExpression.Interp
import HashedExpression.Modeling.Typed
import HashedExpression.Value
import Test.HUnit
import Test.Hspec
import Test.QuickCheck
import Text.Printf
import Var
import Prelude hiding ((**), (^))

sizeReduceFactor :: Int
sizeReduceFactor = 3

shapeToString :: Shape -> String
shapeToString = intercalate "x" . map show

removeDuplicate :: (Ord a) => [a] -> [a]
removeDuplicate = toList . fromList

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

genDouble :: Gen Double
genDouble = arbitrary `suchThat` inSmallRange
  where
    inSmallRange x = x >= 0 && x <= 10

vectorOfDifferent :: Eq a => Int -> Gen a -> Gen [a]
vectorOfDifferent sz gen = foldM f [] [1 .. sz]
  where
    f acc _ = (: acc) <$> gen `suchThat` (not . flip elem acc)

-------------------------------------------------------------------------------
genDimSelector :: Int -> Gen DimSelector
genDimSelector size = do
  let id = elements [0 .. size - 1]
  let step = elements [1 .. size]
  oneof [Range <$> id <*> id <*> step, At <$> id]

genAtSelector :: Int -> Gen DimSelector
genAtSelector size = do
  let id = elements [0 .. size - 1]
  oneof [At <$> id]

genValMapFor :: ExpressionMap -> Gen ValMap
genValMapFor mp = do
  let genVal :: Shape -> Gen Val
      genVal shape = do
        let dbls = vectorOf (product shape) genDouble
        case shape of
          [] -> VScalar <$> genDouble
          [sz] -> V1D . listArray (0, sz - 1) <$> dbls
          [sz1, sz2] -> V2D . listArray ((0, 0), (sz1 - 1, sz2 - 1)) <$> dbls
          [sz1, sz2, sz3] -> V3D . listArray ((0, 0, 0), (sz1 - 1, sz2 - 1, sz3 - 1)) <$> dbls
  vals <-
    (varNodes mp ++ paramNodes mp)
      & mapM
        ( \(name, shape, _) -> do
            val <- genVal shape
            return (name, val)
        )
  return $ Map.fromList vals

primitiveRUntyped :: Shape -> Gen ExprBuilder
primitiveRUntyped shape = do
  dbl <- genDouble
  name <- elements $ map pure ['a' .. 'z']
  let varName = name ++ shapeToString shape
  let parName = "p" ++ name ++ shapeToString shape
  elements
    [ introduceNode (shape, R, Var varName),
      introduceNode (shape, R, Param parName),
      introduceNode (shape, R, Const dbl)
    ]

primitiveCUntyped :: Shape -> Gen ExprBuilder
primitiveCUntyped shape = do
  re <- primitiveRUntyped shape
  im <- primitiveRUntyped shape
  return $ re +: im

genExpUntyped :: Int -> Shape -> ElementType -> Gen ExprBuilder
genExpUntyped qc shape et
  | qc == 0 = case et of
    R -> primitiveRUntyped shape
    C -> primitiveCUntyped shape
  | otherwise =
    let sub = genExpUntyped (qc `div` sizeReduceFactor) shape et
        subScalarR = genExpUntyped (qc `div` sizeReduceFactor) [] R
        subScalarC = genExpUntyped (qc `div` sizeReduceFactor) [] C
        subR = genExpUntyped (qc `div` sizeReduceFactor) shape R
        subC = genExpUntyped (qc `div` sizeReduceFactor) shape C
        subOf shape et = genExpUntyped (qc `div` sizeReduceFactor) shape et
        fromPiecewise = do
          numBranches <- elements [2, 3]
          branches <- vectorOf numBranches sub
          condition <- subR
          marks <- sort <$> vectorOfDifferent (numBranches - 1) arbitrary
          return $ piecewise marks condition branches
        binary op = liftA2 op sub sub
        unary op = op <$> sub
        commonPossibilities =
          [ fromPiecewise,
            binary (+),
            binary (*),
            liftA2 (*.) subScalarR sub,
            unary (negate),
            unary (^ 2)
          ]
        specificElementTypePosibilities
          | et == R =
            [ xRe <$> subC,
              xIm <$> subC
            ]
          | et == C =
            [ liftA2 (+:) subR subR,
              unary ft,
              unary ift,
              liftA2 (*.) subScalarC sub
            ]
        specificShapePossibilities = case shape of
          [] ->
            let fromProjection = do
                  size <- elements [3 .. 9]
                  exp <- subOf [size] et
                  ds <- genAtSelector size
                  return $ project [ds] exp
             in [ do
                    size <- elements [3 .. 9]
                    x <- subOf [size] et
                    y <- subOf [size] et
                    return $ x <.> y,
                  do
                    size1 <- elements [3 .. 9]
                    size2 <- elements [3 .. 9]
                    x <- subOf [size1, size2] et
                    y <- subOf [size1, size2] et
                    return $ x <.> y,
                  fromProjection
                ]
          [n] ->
            let fromRotate = do
                  amount <- elements [- n .. n]
                  exp <- sub
                  return $ rotate [amount] exp
                fromProjectInject = do
                  exp1 <- sub
                  exp2 <- sub
                  ds <- genDimSelector n
                  return $ inject [ds] (project [ds] exp1) exp2
                fromMatrixMul = do
                  m <- elements [3 .. 9]
                  x <- subOf [n, m] et
                  y <- subOf [m] et
                  return $ x ** y
             in [ fromRotate,
                  fromProjectInject,
                  fromMatrixMul
                ]
          [m, n] ->
            let fromRotate = do
                  amount1 <- elements [- m .. m]
                  amount2 <- elements [- n .. n]
                  exp <- sub
                  return $ rotate [amount1, amount2] exp
                fromProjectInject = do
                  exp1 <- sub
                  exp2 <- sub
                  ds1 <- genDimSelector m
                  ds2 <- genDimSelector n
                  return $ inject [ds1, ds2] (project [ds1, ds2] exp1) exp2
                fromMatrixMul = do
                  p <- elements [3 .. 9]
                  x <- subOf [m, p] et
                  y <- subOf [p, n] et
                  return $ x ** y
                fromTranspose = do
                  x <- subOf [n, m] et
                  return $ transpose x
             in [ fromRotate,
                  fromProjectInject,
                  fromMatrixMul,
                  fromTranspose
                ]
     in oneof $ commonPossibilities ++ specificShapePossibilities ++ specificElementTypePosibilities

genExp :: forall d et. (ToShape d, IsElementType et) => Int -> Gen (TypedExpr d et)
genExp size = TypedExpr <$> genExpUntyped size (toShape @d) (toElementType @et)

-------------------------------------------------------------------------------
instance (ToShape d, IsElementType et) => Arbitrary (TypedExpr d et) where
  arbitrary = sized genExp

-------------------------------------------------------------------------------
data Suite d et
  = Suite (TypedExpr d et) ValMap
  deriving (Show)

instance (ToShape d, IsElementType et) => Arbitrary (Suite d et) where
  arbitrary = do
    exp <- arbitrary
    valMap <- genValMapFor (fst $ asRawExpr exp)
    return $ Suite exp valMap

-------------------------------------------------------------------------------
type SuiteScalarR = Suite Scalar R

type SuiteScalarC = Suite Scalar C

type SuiteOneR = Suite '[Default1D] R

type SuiteOneC = Suite '[Default1D] C

type SuiteTwoR = Suite '[Default2D1, Default2D2] R

type SuiteTwoC = Suite '[Default2D1, Default2D2] C

-------------------------------------------------------------------------------
newtype ArbitraryExpr = ArbitraryExpr {unArbitraryExpr :: RawExpr}
  deriving (Show, Ord, Eq)

instance Arbitrary ArbitraryExpr where
  arbitrary = do
    elementType <- elements [R, C]
    shapeLength <- elements [0 .. 2]
    shape <- vectorOf shapeLength $ elements [1 .. 6]
    ArbitraryExpr . buildExpr <$> sized (\sz -> genExpUntyped sz shape elementType)

data XSuite = XSuite RawExpr ValMap deriving (Show, Eq, Ord)

instance Arbitrary XSuite where
  arbitrary = do
    ArbitraryExpr exp <- arbitrary
    valMap <- genValMapFor (fst exp)
    return $ XSuite exp valMap

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

instance (Ix i, Num a) => Num (Array i a) where
  (+) arr1 arr2 = listArray (bounds arr1) $ zipWith (+) (elems arr1) (elems arr2)
  (*) arr1 arr2 = listArray (bounds arr1) $ zipWith (*) (elems arr1) (elems arr2)

instance (Ix i) => ComplexRealOp (Array i Double) (Array i (Complex Double)) where
  (+:) arr1 arr2 = listArray (bounds arr1) $ zipWith (+:) (elems arr1) (elems arr2)
  xRe arr = listArray (bounds arr) $ map xRe (elems arr)
  xIm arr = listArray (bounds arr) $ map xIm (elems arr)
  conjugate arr = listArray (bounds arr) $ map conjugate (elems arr)

-- | Approximable class
class Show a => Approximable a where
  -- | (~=) checks if the two values are within acceptable numerical error
  (~=) :: a -> a -> Bool

  -- | Prettyprint a string produced by show
  prettifyShow :: a -> String

infix 4 ~=

{-
   Calculating the reletive error for two double-precision input arguments.
   It returns the error calculating by |a-b|/ max(|a|.|b|), where a and b are double-precision numbers, so it returns the relative
   error as double.
-}
relativeError :: Double -> Double -> Double
relativeError a b = abs (a - b) / max (abs a) (abs b)

instance Approximable Double where
  (~=) :: Double -> Double -> Bool
  a ~= b
    | abs (a - b) < 1.0e-5 = True
    | a == b = True
    | otherwise = relativeError a b < 0.01
  prettifyShow a
    | abs a < 1e-10 = "0" --  in prettify, inputs less than a small condition value consider as zero,
    | otherwise = printf "%.2f" a --  otherwise, it shows the exact input.

{-
Instance which belongs to approximable class for Complex Double precision inputs.
It takes 2 Complex double precision inputs and apply ~= operation on them.
calculating the error in real and imaginary parts of the complex input seperately.
Returns true if the error is less than the condition value.
-}
instance Approximable (Complex Double) where
  (~=) :: Complex Double -> Complex Double -> Bool -- Using approximable operation ~= for two Complex Double input and returns Bool
  a ~= b = (realPart a ~= realPart b) && (imagPart a ~= imagPart b) --  check for real and imaginary parts seperately
  prettifyShow a =
    prettifyShow (realPart a) ++ " + " ++ prettifyShow (imagPart a) ++ "i" --  Prettyprint a string produced by show

instance Approximable (Array Int Double) where
  (~=) :: Array Int Double -> Array Int Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable (Array Int (Complex Double)) where
  (~=) :: Array Int (Complex Double) -> Array Int (Complex Double) -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable (Array (Int, Int) Double) where
  (~=) :: Array (Int, Int) Double -> Array (Int, Int) Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable (Array (Int, Int) (Complex Double)) where
  (~=) ::
    Array (Int, Int) (Complex Double) ->
    Array (Int, Int) (Complex Double) ->
    Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable (Array (Int, Int, Int) Double) where
  (~=) :: Array (Int, Int, Int) Double -> Array (Int, Int, Int) Double -> Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable (Array (Int, Int, Int) (Complex Double)) where
  (~=) ::
    Array (Int, Int, Int) (Complex Double) ->
    Array (Int, Int, Int) (Complex Double) ->
    Bool
  a ~= b = (bounds a == bounds b) && and (zipWith (~=) (elems a) (elems b))
  prettifyShow a =
    "[" ++ (intercalate ", " . map prettifyShow . elems $ a) ++ "]" --  concatenate "," to seperate the elements of array in prettify

instance Approximable InterpValue where
  VR x ~= VR y = x ~= y
  V1DR x ~= V1DR y = x ~= y
  V2DR x ~= V2DR y = x ~= y
  V3DR x ~= V3DR y = x ~= y
  VC x ~= VC y = x ~= y
  V1DC x ~= V1DC y = x ~= y
  V2DC x ~= V2DC y = x ~= y
  V3DC x ~= V3DC y = x ~= y
  _ ~= _ = False

  prettifyShow a = case a of
    VR x -> prettifyShow x
    V1DR x -> prettifyShow x
    V2DR x -> prettifyShow x
    V3DR x -> prettifyShow x
    VC x -> prettifyShow x
    V1DC x -> prettifyShow x
    V2DC x -> prettifyShow x
    V3DC x -> prettifyShow x

-------------------------------------------------------------------------------
shouldApprox :: (HasCallStack, Approximable a) => a -> a -> Expectation
shouldApprox x y = assertBool msg (x ~= y)
  where
    msg = "Expected: " ++ prettifyShow y ++ "\nGot: " ++ prettifyShow x

infix 1 `shouldApprox`

max_stars_repo_name: McMasterU/HashedExpression

max_issues_repo_name: McMasterU/HashedExpression

max_forks_repo_name: McMasterU/HashedExpression


text:
{-|
Module      : Lib
Description : Main component of your project.
Copyright   : (c) Fabricio Olivetti de Franca, 2021
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com
Stability   : experimental
Portability : POSIX
|-}

module GAPoly
    ( parseFile
    , runGA
    , generateReports
    , Solution(..)
    , Poly(..)
    ) where

import GA
import Fitness
import Dataset
import Report
import Random
import Evolution

import Control.Monad.State.Strict
import System.Random
import Numeric.LinearAlgebra (Vector)
import GHC.TypeLits

gaPoly :: KnownNat n => Double -> Evolution (Poly n)
gaPoly pm = Select generational End
              (Cross crossover
                (Mutate (mutate pm) End))

-- | Run a genetic algorithm with the provided parameters
runGA :: KnownNat n => Int -> Int -> Int -> Int -> Double -> [[Double]] -> Vector Double -> IO ([Double], Solution (Poly n))
runGA nTerms nVars it nPop pm xss ys =
  do g <- newStdGen
     let fitness   = evalFitness nTerms xss ys
         createSol = createRndSolution (nVars*nTerms)
         gens      = runEvolution createSol fitness it nPop (gaPoly pm)
     evalStateT gens g

max_stars_repo_name: folivetti/gapoly

max_issues_repo_name: folivetti/gapoly

max_forks_repo_name: folivetti/gapoly


text:
{-# LANGUAGE BangPatterns, FlexibleContexts #-}
-- |
-- Module    : Statistics.Transform
-- Copyright : (c) 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Fourier-related transformations of mathematical functions.
--
-- These functions are written for simplicity and correctness, not
-- speed.  If you need a fast FFT implementation for your application,
-- you should strongly consider using a library of FFTW bindings
-- instead.

module Statistics.Transform
    (
    -- * Type synonyms
      CD
    -- * Discrete cosine transform
    , dct
    , dct_
    , idct
    , idct_
    -- * Fast Fourier transform
    , fft
    , ifft
    ) where

import Control.Monad (when)
import Control.Monad.ST (ST)
import Data.Bits (shiftL, shiftR)
import Data.Complex (Complex(..), conjugate, realPart)
import Numeric.SpecFunctions (log2)
import qualified Data.Vector.Generic         as G
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Unboxed         as U
import qualified Data.Vector                 as V

type CD = Complex Double

-- | Discrete cosine transform (DCT-II).
dct :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v Double -> v Double
dct = dctWorker . G.map (:+0)
{-# INLINABLE  dct #-}
{-# SPECIAlIZE dct :: U.Vector Double -> U.Vector Double #-}
{-# SPECIAlIZE dct :: V.Vector Double -> V.Vector Double #-}

-- | Discrete cosine transform (DCT-II). Only real part of vector is
--   transformed, imaginary part is ignored.
dct_ :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v CD -> v Double
dct_ = dctWorker . G.map (\(i :+ _) -> i :+ 0)
{-# INLINABLE  dct_ #-}
{-# SPECIAlIZE dct_ :: U.Vector CD -> U.Vector Double #-}
{-# SPECIAlIZE dct_ :: V.Vector CD -> V.Vector Double#-}

dctWorker :: (G.Vector v CD, G.Vector v Double, G.Vector v Int) => v CD -> v Double
{-# INLINE dctWorker #-}
dctWorker xs
  -- length 1 is special cased because shuffle algorithms fail for it.
  | G.length xs == 1 = G.map ((2*) . realPart) xs
  | vectorOK xs      = G.map realPart $ G.zipWith (*) weights (fft interleaved)
  | otherwise        = error "Statistics.Transform.dct: bad vector length"
  where
    interleaved = G.backpermute xs $ G.enumFromThenTo 0 2 (len-2) G.++
                                     G.enumFromThenTo (len-1) (len-3) 1
    weights = G.cons 2 . G.generate (len-1) $ \x ->
              2 * exp ((0:+(-1))*fi (x+1)*pi/(2*n))
      where n = fi len
    len = G.length xs

-- | Inverse discrete cosine transform (DCT-III). It's inverse of
-- 'dct' only up to scale parameter:
--
-- > (idct . dct) x = (* length x)
idct :: (G.Vector v CD, G.Vector v Double) => v Double -> v Double
idct = idctWorker . G.map (:+0)
{-# INLINABLE  idct #-}
{-# SPECIAlIZE idct :: U.Vector Double -> U.Vector Double #-}
{-# SPECIAlIZE idct :: V.Vector Double -> V.Vector Double #-}

-- | Inverse discrete cosine transform (DCT-III). Only real part of vector is
--   transformed, imaginary part is ignored.
idct_ :: (G.Vector v CD, G.Vector v Double) => v CD -> v Double
idct_ = idctWorker . G.map (\(i :+ _) -> i :+ 0)
{-# INLINABLE  idct_ #-}
{-# SPECIAlIZE idct_ :: U.Vector CD -> U.Vector Double #-}
{-# SPECIAlIZE idct_ :: V.Vector CD -> V.Vector Double #-}

idctWorker :: (G.Vector v CD, G.Vector v Double) => v CD -> v Double
{-# INLINE idctWorker #-}
idctWorker xs
  | vectorOK xs = G.generate len interleave
  | otherwise   = error "Statistics.Transform.dct: bad vector length"
  where
    interleave z | even z    = vals `G.unsafeIndex` halve z
                 | otherwise = vals `G.unsafeIndex` (len - halve z - 1)
    vals = G.map realPart . ifft $ G.zipWith (*) weights xs
    weights
      = G.cons n
      $ G.generate (len - 1) $ \x -> 2 * n * exp ((0:+1) * fi (x+1) * pi/(2*n))
      where n = fi len
    len = G.length xs

-- | Inverse fast Fourier transform.
ifft :: G.Vector v CD => v CD -> v CD
ifft xs
  | vectorOK xs = G.map ((/fi (G.length xs)) . conjugate) . fft . G.map conjugate $ xs
  | otherwise   = error "Statistics.Transform.ifft: bad vector length"
{-# INLINABLE  ifft #-}
{-# SPECIAlIZE ifft :: U.Vector CD -> U.Vector CD #-}
{-# SPECIAlIZE ifft :: V.Vector CD -> V.Vector CD #-}

-- | Radix-2 decimation-in-time fast Fourier transform.
fft :: G.Vector v CD => v CD -> v CD
fft v | vectorOK v  = G.create $ do mv <- G.thaw v
                                    mfft mv
                                    return mv
      | otherwise   = error "Statistics.Transform.fft: bad vector length"
{-# INLINABLE  fft #-}
{-# SPECIAlIZE fft :: U.Vector CD -> U.Vector CD #-}
{-# SPECIAlIZE fft :: V.Vector CD -> V.Vector CD #-}

-- Vector length must be power of two. It's not checked
mfft :: (M.MVector v CD) => v s CD -> ST s ()
{-# INLINE mfft #-}
mfft vec = bitReverse 0 0
 where
  bitReverse i j | i == len-1 = stage 0 1
                 | otherwise  = do
    when (i < j) $ M.swap vec i j
    let inner k l | k <= l    = inner (k `shiftR` 1) (l-k)
                  | otherwise = bitReverse (i+1) (l+k)
    inner (len `shiftR` 1) j
  stage l !l1 | l == m    = return ()
              | otherwise = do
    let !l2 = l1 `shiftL` 1
        !e  = -6.283185307179586/fromIntegral l2
        flight j !a | j == l1   = stage (l+1) l2
                    | otherwise = do
          let butterfly i | i >= len  = flight (j+1) (a+e)
                          | otherwise = do
                let i1 = i + l1
                xi1 :+ yi1 <- M.read vec i1
                let !c = cos a
                    !s = sin a
                    d  = (c*xi1 - s*yi1) :+ (s*xi1 + c*yi1)
                ci <- M.read vec i
                M.write vec i1 (ci - d)
                M.write vec i (ci + d)
                butterfly (i+l2)
          butterfly j
    flight 0 0
  len = M.length vec
  m   = log2 len

----------------------------------------------------------------
-- Helpers
----------------------------------------------------------------

fi :: Int -> CD
fi = fromIntegral

halve :: Int -> Int
halve = (`shiftR` 1)

vectorOK :: G.Vector v a => v a -> Bool
{-# INLINE vectorOK #-}
vectorOK v = (1 `shiftL` log2 n) == n where n = G.length v

max_stars_repo_name: StefanHubner/statistics

max_issues_repo_name: StefanHubner/statistics

max_forks_repo_name: StefanHubner/statistics


text:
import Data.Complex
import Data.Array
import Data.Ratio
import qualified Data.Map as M

fft :: [Complex Double] -> [Complex Double]
fft x = let n = length x
            i = 0 :+ 1
            w = M.fromList [(k%n, exp ((-2)*pi*i*(fromIntegral k)/(fromIntegral n)) ) | k<-[0..n-1]]
            arr = fft' n w (listArray (0,n-1) x)
        in [arr!k | k<-[0..n-1]]
  where
  fft' 1 _ x = x
  fft' n w x = let n2 = div n 2
                   e = fft' n2 w (listArray (0, n2-1) [x!k | k<-[0,2..n-1]])
                   o = fft' n2 w (listArray (0, n2-1) [x!k | k<-[1,3..n-1]])
               in array (0, n-1) $ concat [[(k, e!k + o!k * w M.!(k%n)),
                                            (k + n2, e!k - o!k * w M.!(k%n))]
                                                               | k <- [0..n2-1]]

main = do
  print $ fft [0,1,2,3]

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
module Main where

import Codec.Picture
import Codec.Picture.Gif
import Codec.Picture.Png
import Codec.Picture.Types
import Control.Parallel
import Control.Parallel.Strategies
import Data.Complex
import Data.List.Split
import Data.Vector.Storable
import System.Environment

main = do
  putStrLn "Resolution: "
  size <- read <$> getLine
  putStrLn "Iterations: "
  iters <- read <$> getLine
  writePng "out.png" $ genFractal size iters 2

genFractal :: Int -> Int -> Complex Double -> Image Pixel8
genFractal size iters comp = Image size size $ fromList pixels
  where
    pixels =
      [genPixel comp (i `mod` size) (i `div` size) | i <- [0 .. size * size `div` 2]] `using`
      parList rdeepseq :: [Pixel8]
    genPixel comp x y =
      fromIntegral $ div 0xff iters * (iters - boundedUntil comp) :: Pixel8
      where
        boundedUntil =
          isBounded iters $
          fromIntegral (x - div size 2) / (fromIntegral size / 3) :+
          fromIntegral (y - div size 2) / (fromIntegral size / 3)

isBounded iters c comp = isBounded' 1 c
  where
    isBounded' i z
      | i == iters = i
      | magnitude z > 2 = i
      | otherwise = isBounded' (i + 1) $ (z ** comp) + c

max_stars_repo_name: IanPoggers/fractal

max_issues_repo_name: IanPoggers/fractal

max_forks_repo_name: IanPoggers/fractal


text:
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE RecordWildCards       #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE FlexibleContexts      #-}
{-|
Module      : Grenade.Layers.Deconvolution
Description : Deconvolution layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental

A deconvolution layer is in many ways a convolution layer in reverse.
It learns a kernel to apply to each pixel location, spreading it out
into a larger layer.

This layer is important for image generation tasks, such as GANs on
images.
-}
module Grenade.Layers.Deconvolution (
    Deconvolution (..)
  , Deconvolution' (..)
  ) where

import           Data.Maybe
import           Data.Proxy
import           Data.Serialize
import           Data.Singletons.TypeLits hiding (natVal)

import           GHC.TypeLits

import Control.DeepSeq (NFData (..))

import           Numeric.LinearAlgebra hiding ( uniformSample, konst )
import qualified Numeric.LinearAlgebra as LA
import           Numeric.LinearAlgebra.Static hiding ((|||), build, toRows)

import           Grenade.Core
import           Grenade.Layers.Internal.Convolution
import           Grenade.Layers.Internal.Update

-- | A Deconvolution layer for a neural network.
--   This uses the im2col Convolution trick popularised by Caffe.
--
--   The Deconvolution layer is a way of spreading out a single response
--   into a larger image, and is useful in generating images.
--
data Deconvolution :: Nat -- Number of channels, for the first layer this could be RGB for instance.
                   -> Nat -- Number of filters, this is the number of channels output by the layer.
                   -> Nat -- The number of rows in the kernel filter
                   -> Nat -- The number of column in the kernel filter
                   -> Nat -- The row stride of the Deconvolution filter
                   -> Nat -- The columns stride of the Deconvolution filter
                   -> * where
  Deconvolution :: ( KnownNat channels
                   , KnownNat filters
                   , KnownNat kernelRows
                   , KnownNat kernelColumns
                   , KnownNat strideRows
                   , KnownNat strideColumns
                   , KnownNat kernelFlattened
                   , kernelFlattened ~ (kernelRows * kernelColumns * filters))
                 => !(L kernelFlattened channels) -- The kernel filter weights
                 -> !(L kernelFlattened channels) -- The last kernel update (or momentum)
                 -> Deconvolution channels filters kernelRows kernelColumns strideRows strideColumns

instance NFData (Deconvolution c f k k' s s') where
  rnf (Deconvolution a b) = rnf a `seq` rnf b `seq` ()

data Deconvolution' :: Nat -- Number of channels, for the first layer this could be RGB for instance.
                    -> Nat -- Number of filters, this is the number of channels output by the layer.
                    -> Nat -- The number of rows in the kernel filter
                    -> Nat -- The number of column in the kernel filter
                    -> Nat -- The row stride of the Deconvolution filter
                    -> Nat -- The columns stride of the Deconvolution filter
                    -> * where
  Deconvolution' :: ( KnownNat channels
                  , KnownNat filters
                  , KnownNat kernelRows
                  , KnownNat kernelColumns
                  , KnownNat strideRows
                  , KnownNat strideColumns
                  , KnownNat kernelFlattened
                  , kernelFlattened ~ (kernelRows * kernelColumns * filters))
               => !(L kernelFlattened channels) -- The kernel filter gradient
               -> Deconvolution' channels filters kernelRows kernelColumns strideRows strideColumns

instance NFData (Deconvolution' c f k k' s s') where
  rnf (Deconvolution' a) = rnf a `seq` ()

instance Show (Deconvolution c f k k' s s') where
  show (Deconvolution a _) = renderConv a
    where
      renderConv mm =
        let m  = extract mm
            ky = fromIntegral $ natVal (Proxy :: Proxy k)
            rs = LA.toColumns m
            ms = map (take ky) $ toLists . reshape ky <$> rs

            render n'  | n' <= 0.2  = ' '
                       | n' <= 0.4  = '.'
                       | n' <= 0.6  = '-'
                       | n' <= 0.8  = '='
                       | otherwise =  '#'

            px = (fmap . fmap . fmap) render ms
        in unlines $ foldl1 (zipWith (\a' b' -> a' ++ "   |   " ++ b')) $ px

instance ( KnownNat c
         , KnownNat f
         , KnownNat k
         , KnownNat k'
         , KnownNat s
         , KnownNat s'
         , KnownNat ((k * k') * f)
         , KnownNat ((k * k') * c)
         , KnownNat (c * ((k * k') * f))) => RandomLayer (Deconvolution c f k k' s s') where
  createRandomWith m gen = do
    wN <- getRandomMatrix i i m gen
    let mm = konst 0
    return $ Deconvolution wN mm
    where i = natVal (Proxy :: Proxy ((k * k') * c))

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * filters)
         ) => UpdateLayer (Deconvolution channels filters kernelRows kernelColumns strideRows strideColumns) where
  type Gradient (Deconvolution channels filters kernelRows kernelColumns strideRows strideColumns) = (Deconvolution' channels filters kernelRows kernelColumns strideRows strideColumns)
  runUpdate LearningParameters {..} (Deconvolution oldKernel oldMomentum) (Deconvolution' kernelGradient) =
    let (newKernel, newMomentum) = descendMatrix learningRate learningMomentum learningRegulariser oldKernel kernelGradient oldMomentum
    in Deconvolution newKernel newMomentum

instance ( KnownNat channels
         , KnownNat filters
         , KnownNat kernelRows
         , KnownNat kernelColumns
         , KnownNat strideRows
         , KnownNat strideColumns
         , KnownNat (kernelRows * kernelColumns * filters)
         ) => Serialize (Deconvolution channels filters kernelRows kernelColumns strideRows strideColumns) where
  put (Deconvolution w _) = putListOf put . toList . flatten . extract $ w
  get = do
      let f  = fromIntegral $ natVal (Proxy :: Proxy channels)
      wN    <- maybe (fail "Vector of incorrect size") return . create . reshape f . LA.fromList =<< getListOf get
      let mm = konst 0
      return $ Deconvolution wN mm

-- | A two dimentional image may have a Deconvolution filter applied to it
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , ((inputRows - 1) * strideRows) ~ (outputRows - kernelRows)
         , ((inputCols - 1) * strideCols) ~ (outputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * filters)
         , KnownNat (outputRows * filters)
         ) => Layer (Deconvolution 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) where
  type Tape (Deconvolution 1 filters kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    runForwards c (S3D input :: S ('D3 inputRows inputCols 1))

  runBackwards c tape grads =
    case runBackwards c tape grads of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

-- | A two dimentional image may have a Deconvolution filter applied to it
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , ((inputRows - 1) * strideRows) ~ (outputRows - kernelRows)
         , ((inputCols - 1) * strideCols) ~ (outputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * 1)
         , KnownNat (outputRows * 1)
         ) => Layer (Deconvolution 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) where
  type Tape (Deconvolution 1 1 kernelRows kernelCols strideRows strideCols) ('D2 inputRows inputCols) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols 1)
  runForwards c (S2D input) =
    case runForwards c (S3D input :: S ('D3 inputRows inputCols 1)) of
      (tps, S3D fore :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D fore)

  runBackwards c tape (S2D grads) =
    case runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1)) of
      (c', S3D back :: S ('D3 inputRows inputCols 1)) ->  (c', S2D back)

-- | A two dimentional image may have a Deconvolution filter applied to it
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , ((inputRows - 1) * strideRows) ~ (outputRows - kernelRows)
         , ((inputCols - 1) * strideCols) ~ (outputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * 1)
         , KnownNat (outputRows * 1)
         , KnownNat channels
         ) => Layer (Deconvolution channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) where
  type Tape (Deconvolution channels 1 kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D2 outputRows outputCols) = S ('D3 inputRows inputCols channels)
  runForwards c input =
    case runForwards c input of
      (tps, S3D fore :: S ('D3 outputRows outputCols 1)) ->  (tps, S2D fore)

  runBackwards c tape (S2D grads) =
    runBackwards c tape (S3D grads :: S ('D3 outputRows outputCols 1))

-- | A three dimensional image (or 2d with many channels) can have
--   an appropriately sized Deconvolution filter run across it.
instance ( KnownNat kernelRows
         , KnownNat kernelCols
         , KnownNat filters
         , KnownNat strideRows
         , KnownNat strideCols
         , KnownNat inputRows
         , KnownNat inputCols
         , KnownNat outputRows
         , KnownNat outputCols
         , KnownNat channels
         , ((inputRows - 1) * strideRows) ~ (outputRows - kernelRows)
         , ((inputCols - 1) * strideCols) ~ (outputCols - kernelCols)
         , KnownNat (kernelRows * kernelCols * filters)
         , KnownNat (outputRows * filters)
         ) => Layer (Deconvolution channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) where

  type Tape (Deconvolution channels filters kernelRows kernelCols strideRows strideCols) ('D3 inputRows inputCols channels) ('D3 outputRows outputCols filters) = S ('D3 inputRows inputCols channels)

  runForwards (Deconvolution kernel _) (S3D input) =
    let ex = extract input
        ek = extract kernel
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)

        c  = vid2col 1 1 1 1 ix iy ex

        mt = c LA.<> tr ek

        r  = col2vid kx ky sx sy ox oy mt
        rs = fromJust . create $ r
    in  (S3D input, S3D rs)
  runBackwards (Deconvolution kernel _) (S3D input) (S3D dEdy) =
    let ex = extract input
        ix = fromIntegral $ natVal (Proxy :: Proxy inputRows)
        iy = fromIntegral $ natVal (Proxy :: Proxy inputCols)
        kx = fromIntegral $ natVal (Proxy :: Proxy kernelRows)
        ky = fromIntegral $ natVal (Proxy :: Proxy kernelCols)
        sx = fromIntegral $ natVal (Proxy :: Proxy strideRows)
        sy = fromIntegral $ natVal (Proxy :: Proxy strideCols)
        ox = fromIntegral $ natVal (Proxy :: Proxy outputRows)
        oy = fromIntegral $ natVal (Proxy :: Proxy outputCols)

        c  = vid2col 1 1 1 1 ix iy ex

        eo = extract dEdy
        ek = extract kernel

        vs = vid2col kx ky sx sy ox oy eo

        kN = fromJust . create . tr $ tr c LA.<> vs

        dW = vs LA.<> ek

        xW = col2vid 1 1 1 1 ix iy dW
    in  (Deconvolution' kN, S3D . fromJust . create $ xW)

-------------------- GNum instances --------------------

instance (KnownNat strideCols,KnownNat strideRows,KnownNat kernelCols,KnownNat kernelRows,KnownNat filters,KnownNat channels,KnownNat ((kernelRows * kernelCols) * filters),KnownNat
                          ((kernelRows * kernelCols) * channels)) => GNum (Deconvolution channels filters kernelRows kernelCols strideRows strideCols) where
  n |* (Deconvolution w m) = Deconvolution (fromRational n * w) m
  (Deconvolution w m) |+ (Deconvolution w2 m2)  = Deconvolution (fromRational 0.5 * (w+w2)) (fromRational 0.5 * (m+m2))
  gFromRational r = Deconvolution (fromRational r) (fromRational r)

instance (KnownNat strideCols,KnownNat strideRows,KnownNat kernelCols,KnownNat kernelRows,KnownNat filters,KnownNat channels,KnownNat ((kernelRows * kernelCols) * filters),KnownNat
                          ((kernelRows * kernelCols) * channels)) => GNum (Deconvolution' channels filters kernelRows kernelCols strideRows strideCols) where
  _ |* (Deconvolution' g) = Deconvolution' g
  (Deconvolution' g) |+ (Deconvolution' g2)  = Deconvolution' (fromRational 0.5 * (g+g2))
  gFromRational r = Deconvolution' (fromRational r)

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


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

module Lib
    ( trilinear
    , resample
    , affine
    , array_correlation
    , fisher_transform
    , naive_12p_registration
    ) where

import Data.Array.Repa
import qualified Data.Array.Repa as Repa
import Data.Array.Repa.Algorithms.Matrix
import Statistics.Sample
import Data.Array.Repa.Repr.Unboxed

type DIM3D = (Z :. Double :. Double :. Double)

trilinear :: (Unbox a, Floating a, Fractional a) => Array U DIM3 a -> DIM3D -> a
trilinear !arr !(Z :. x :. y :. z) =
  let
    x0 = floor x
    x1 = ceiling x :: Int
    y0 = floor y
    y1 = ceiling y :: Int
    z0 = floor z
    z1 = ceiling z :: Int
    in let
      v000 = arr ! (Z :. x0 :. y0 :. z0)
      v100 = arr ! (Z :. x1 :. y0 :. z0)
      v010 = arr ! (Z :. x0 :. y1 :. z0)
      v110 = arr ! (Z :. x1 :. y1 :. z0)
      v001 = arr ! (Z :. x0 :. y0 :. z1)
      v101 = arr ! (Z :. x1 :. y0 :. y1)
      v011 = arr ! (Z :. x0 :. y1 :. y1)
      v111 = arr ! (Z :. x1 :. y1 :. y1)
      xd = if x0 /= x1 then
        ((fromIntegral x0) - (realToFrac x)) / (fromIntegral (x0 - x1)) else
        0
      yd = if y0 /= y1 then
        ((fromIntegral y0) - (realToFrac y)) / (fromIntegral (y0 - y1)) else
        0
      zd = if z0 /= z1 then
        ((fromIntegral z0) - realToFrac z) / (fromIntegral (z0 - z1)) else
        0
      in let
        c00 = v000 * (1 - xd) + v100 * xd
        c01 = v001 * (1 - xd) + v101 * xd
        c10 = v010 * (1 - xd) + v110 * xd
        c11 = v011 * (1 - xd) + v111 * xd
        in let
          c0 = c00 * (1 - yd) + c10 * yd
          c1 = c01 * (1 - yd) + c11 * yd
          in  c0 * (1 - zd) + c1 * zd

resample ::
  (Floating a, Floating b, Unbox a, Unbox b, Fractional a, Fractional b) =>
  Array U DIM3 a ->
  (DIM3 -> DIM3D) -> --fun from coordinate system to coordinate system
  (Array U DIM3 a -> DIM3D -> b) -> --interpolator
  Array D DIM3 b

resample !arr !xfm !interp =
  Repa.traverse arr id xfm_and_interp
  where
    xfm_and_interp = \_ d -> ((interp arr (xfm d)))

affine ::
  (Unbox a, Floating a, Fractional a, Real a) =>
  DIM3 -> Array U DIM2 a -> DIM3D

affine !(Z :. x :. y  :. z) !xfm =
  (Z :. x' :. y' :. z')
  where
    [x',y',z',_] = (toList (mmultS (computeS $ Repa.map realToFrac xfm) pnt_array))
      where
        pnt_array = (fromListUnboxed
                     (Z :. (3 :: Int) :. (1 :: Int))
                     (fmap fromIntegral [x,y,z])) -- :: Array U DIM2 Double

array_correlation ::
  (Unbox a, Real a, Fractional a) => Array U DIM3 a -> Array U DIM3 a -> a
array_correlation !arr1 !arr2 =
  realToFrac
   (correlation
    (toUnboxed
     (computeS
      (Repa.zipWith coerce_and_tuple
         arr1 arr2))))
  where
     coerce_and_tuple x y = (realToFrac x, realToFrac y)

fisher_transform = atanh

-- Test code
naive_12p_registration ::
  (Unbox a, Floating a, Fractional a, Real a) =>
  Array U DIM3 a -> [a] -> a

naive_12p_registration !arr !params =
  array_correlation arr (computeS $ resample arr (flip affine xfm) trilinear)
  where
    xfm = fromListUnboxed (Z :. 4 :. 4) params


text:
module Trade.Algorithm.MovingAverage where

import Statistics.Sample (mean)

import qualified Data.Vector as Vec
import Data.Vector (Vector)

import Data.List (scanl')

newtype WindowSize = WindowSize Int deriving (Show, Eq, Ord)

movingAvgL :: WindowSize -> [Double] -> [Double]
movingAvgL (WindowSize k) lst =
  let (h, t) = splitAt k lst
  in map (/ fromIntegral k) $ scanl' (+) (sum h) $ zipWith (-) t lst

mavgL :: WindowSize -> [Double] -> [(Double, Double)]
mavgL w@(WindowSize k) xs =
  let j = fromIntegral k -- / 2
  in zip [j, j+1 ..] (movingAvgL w xs)

movingAvgV :: WindowSize -> Vector Double -> Vector Double
movingAvgV (WindowSize k) v =
  let len = Vec.length v
      f a = Vec.slice a k v
      idx = Vec.generate (len-k) id
  in Vec.map (mean . f) idx

mavgV :: WindowSize -> Vector Double -> Vector (Double, Double)
mavgV (WindowSize k) v =
  let len = Vec.length v
      f a = (fromIntegral (k+a), mean (Vec.slice a k v))
  in Vec.generate (len-k) f

mavg3 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg3 = mavg (WindowSize 3)

mavg5 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg5 = mavg (WindowSize 5)

mavg8 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg8 = mavg (WindowSize 8)

mavg12 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg12 = mavg (WindowSize 12)

mavg15 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg15 = mavg (WindowSize 15)

mavg20 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg20 = mavg (WindowSize 20)

mavg30 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg30 = mavg (WindowSize 30)

mavg50 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg50 = mavg (WindowSize 50)

mavg100 :: MovingAverage vec => vec Double -> vec (Double, Double)
mavg100 = mavg (WindowSize 100)

mavgBar :: WindowSize -> Vector (t, Double) -> Vector (t, Double)
mavgBar w@(WindowSize k) vs =
  let (ts, ys) = Vec.unzip vs
      ts' = Vec.drop k ts
      ys' = movingAverage w ys
  in Vec.zip ts' ys'

class MovingAverage vec where
  movingAverage :: WindowSize -> vec Double -> vec Double
  mavg :: WindowSize -> vec Double -> vec (Double, Double)

instance MovingAverage [] where
  movingAverage = movingAvgL
  mavg = mavgL

instance MovingAverage Vector where
  movingAverage = movingAvgV
  mavg = mavgV

max_stars_repo_name: fphh/trade

max_issues_repo_name: fphh/trade

max_forks_repo_name: fphh/trade


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

module SMC where

import PPL ( Sample(sample), Infer(score, observe, logscore), Model, runModel, sampleWeighted, WeightedT )
import Randomness ( Random, Weighted, Weight, parallel )
import Control.Monad ( ap, replicateM, zipWithM )
import Data.Either ( partitionEithers )
import Particle ( ParticleT(..), reweight )
import Control.Monad.Trans.Class (lift)
import Data.Monoid (Product(..))
import Control.Monad.Trans.Writer (WriterT, tell, runWriterT)
import Numeric.Log (Log(..))
import Distribution ( Dist (samplefrom), logpdf, Uniform (Uniform) )
import Statistics.Resampling (Bootstrap(resamples))
import Control.Monad.ST (ST)
import Control.Parallel (pseq)
import Data.List ( partition, sortOn, groupBy )
import Debug.Trace ( trace )
import Control.Monad.RWS (gets)

type Address = Int

instance (Monoid w, Sample m) => Sample (ParticleT l w m) where
  sample = Compute . fmap return . sample

instance Sample m => Infer (ParticleT l (Product Weight) m) where
  score w = Done (Product $ Exp $ log w) ()
  logscore w = Done (Product w) ()
  observe d x = Done (Product $ logpdf d x) ()

-- SMC with seperate score and resample statements
type SMCModel = ParticleT () (Product Weight) Random
type AddressSMC = ParticleT Address (Product Weight) Random

avg :: [Weighted Double] -> Double
avg xs = let total = sum (map snd xs) in
         exp $ ln $ sum [Exp (log a) * w / total | (a, w) <- xs]

stats :: Functor f => f [Weighted (Double, Double) ] -> f Double
stats = fmap (avg . map (\((a,b),c) -> (a,c)))

weight :: ParticleT l (Product Weight) m a -> Weight
weight (Done w _) = getProduct w
weight (Resample l w _) = getProduct w
weight (Compute _) = error "Compute has no weight"

-- A stochastic monad is one that can perform some sort of stochastic compute step
class Monad m => Stochastic m  where
  compute :: m a -> Random (m a)

instance Stochastic (ParticleT l (Product Weight) Random) where
  -- run a computation until it reaches a resample or is done
  compute (Compute ma) = ma >>= compute
  compute p = return p

class Stochastic m => Resampleable m where
  resample :: m ()

instance Resampleable SMCModel where
  resample = Resample () mempty $ return ()

class Infer m => SMC m where
  resampleSMC :: Int -> [m a] -> Random [Weighted a]
  normalisingConstant :: Int -> [m a] -> Random Weight

class Infer m => SampleWeighted m where
  doResample :: [m a] -> Random [m a]

-- SMC where every resample is explicit
newtype ExplicitResample a = ExplicitResample { unexplicit :: SMCModel a } deriving (Functor, Applicative, Monad, Sample, Resampleable, Infer, SMC)

-- SMC where every score triggers a resample
newtype AlwaysResample a = AlwaysResample { unalways :: SMCModel a } deriving (Functor, Applicative, Monad, Sample, Resampleable, SMC)

-- SMC where weight zero particles are resampled
newtype AliveFilter a = AliveFilter { unalive :: SMCModel a } deriving (Functor, Applicative, Monad, Sample, Resampleable, SMC)

-- SMC where subsets of particles are resampled based on addresses
newtype AlphaSMC a = AlphaSMC { unalpha :: AddressSMC a } deriving (Functor, Applicative, Monad, Sample, Infer)

-- SMC where subsets of particles are resampled based on addresses
newtype Adaptive a = Adaptive { unadapt :: SMCModel a } deriving (Functor, Applicative, Monad, Sample, Resampleable, Show)

-- SMC with Startified Resampling
newtype Stratified a = Stratified { unstrat :: SMCModel a } deriving (Functor, Applicative, Monad, Sample, Resampleable)

-- SMC with Subset Resampling
newtype Labelled a = Labelled { unlab :: AddressSMC a } deriving (Functor, Applicative, Monad, Sample, Infer)

-- these all have implicit resamples:
instance Infer AliveFilter where
  score w = AliveFilter (score w >> resample)
  observe d x = AliveFilter (observe d x >> resample)

instance Infer Adaptive where
  score w = Adaptive (score w >> resample)
  observe d x = Adaptive (observe d x >> resample)

instance Infer AlwaysResample where
  score w = AlwaysResample (score w >> resample)
  observe d x = AlwaysResample (observe d x >> resample)

instance Infer Stratified where
  score w = Stratified (score w >> resample)
  observe d x = Stratified (observe d x >> resample)

-- Generalized Newtype Deriving isn't quite enough to do these
instance Stochastic ExplicitResample where
  compute (ExplicitResample a) = ExplicitResample <$> compute a

instance Stochastic Labelled where
  compute (Labelled a) = Labelled <$> compute a

instance Stochastic AlwaysResample where
  compute (AlwaysResample a) = AlwaysResample <$> compute a

instance Stochastic Adaptive where
  compute (Adaptive a) = Adaptive <$> compute a

instance Stochastic Stratified where
  compute (Stratified a) = Stratified <$> compute a

resampleLab :: [Labelled a] -> Random [Labelled a]
resampleLab ms = do
    -- do all the stochastic compuations so everything
    -- is either done or at a resample statement
    ms' <- parallel $ map (compute . unlab) ms
    -- check if all particles are done...
    if all isDone ms'
      -- ... in which case return
      then return $ map Labelled ms'
      -- otherwise we do resampling
      else do
        -- we need to prune the resample we dealt with and reset weights
        resampled <- labResample $ map getweight ms'
        -- recursive call to continue inference
        resampleLab $ map Labelled resampled

labResample :: [Weighted (AddressSMC a)] -> Random [AddressSMC a]
labResample ms = let grouped = groupBy (\(x,_) (y,_) -> address x == address y) ms in
    fmap concat $ parallel $ map stratResample grouped

resampleStrat :: [Stratified a] -> Random [Stratified a]
resampleStrat ms = do
    -- do all the stochastic compuations so everything
    -- is either done or at a resample statement
    ms' <- parallel $ map (compute . unstrat) ms
    -- check if all particles are done...
    if all isDone ms'
      -- ... in which case return
      then return $ map Stratified ms'
      -- otherwise we do resampling
      else do
        -- we need to prune the resample we dealt with and reset weights
        resampled <- stratResample $ map getweight ms'
        -- recursive call to continue inference
        resampleStrat $ map Stratified resampled

getSamples :: [Weight] -> [Weighted a] -> [a]
getSamples [] _ = []
getSamples (r:rs) values = let (y:ys) = dropWhile (\(x,w) -> w < r) values
                           in fst y : getSamples rs (y:ys)

residual :: [Weighted a] -> ([a], [Weighted a])
residual [] = ([], [])
residual ((a, w) : rest) = let (xs, wxs) = residual rest
                           in let w' = floor w
                           in (replicate w' a ++ xs, (a, w - fromIntegral w') : wxs)

stratResample :: Functor m => [Weighted (ParticleT l (Product Weight) m a)] -> Random [ParticleT l (Product Weight) m a]
stratResample xs = do
    -- get sum of weights
    let weight = sum $ map snd xs
    let n = length xs
    let avgweight = weight / Exp (log $ fromIntegral n)
    -- residuals
    let (guaranteed, rest) = residual [(a, w / avgweight) | (a,w) <- xs]
    let n' = n - length guaranteed
    -- rest
    let weighted = let (as, ws) = unzip rest in zip as (scanl1 (+) ws)
    let avgrest = snd (last weighted) / Exp (log $ fromIntegral n')

    -- list of n uniform[i-1/n,i/n] random values
    let between a b = do
          r <- samplefrom (Uniform (fromIntegral a) (fromIntegral b))
          return $ avgrest * Exp (log r)

    rs <- zipWithM between [0..n'-1] [1..n']
    return $ map (reweight $ Product avgweight) (getSamples rs weighted ++ guaranteed)
{-
stratResample :: [Weighted (SMCModel a)] -> Random [SMCModel a]
stratResample xs = do
    -- [(a1, w1), (a2, w2), ...] becomes [(a1, w1), (a2, w1+w2), ...]
    let weighted = let (as, ws) = unzip xs in zip as (scanl1 (+) ws)
    -- get sum of weights
    let weight = snd $ last weighted
    let n = length weighted
    let avgweight = weight / Exp (log $ fromIntegral n)

    -- list of n uniform[i-1/n,i/n] random values
    let between a b = do
          r <- samplefrom (Uniform (fromIntegral a) (fromIntegral b))
          return $ avgweight * Exp (log r)
    rs <- zipWithM between [0..n-1] [1..n]

    -- list of k weighted samples
    let result = getSamples rs weighted
    return $ map (reweight $ Product avgweight) result
-}

class Infer m => Addressable m where
  resample' :: Address -> m ()

  score' :: Address -> Double -> m ()
  score' a w = score w >> resample' a

  observe' :: Dist d v => Address -> d -> v -> m ()
  observe' a d v = observe d v >> resample' a

-- subset SMC specific functions
instance Addressable Labelled where
  resample' a = Labelled $ Resample a mempty $ return ()

-- alpha-SMC specific functions
instance Addressable AlphaSMC where
  resample' a = AlphaSMC $ Resample a mempty $ return ()

-- get address from a Dome or Resample
address :: Bounded l => ParticleT l w m a -> l
address (Resample a _ _) = a
address (Done _ _) = maxBound
address (Compute _ ) = error "compute has no address"

-- helper function used by alpha-SMC
resampleMinSubset :: [AddressSMC a] -> Random [AddressSMC a]
resampleMinSubset ps = do
  let mini = minimum $ map address ps
  let (minima, rest) = partition ((== mini) . address) ps
  let w = sum $ map weight minima
  let n = length minima
  minima' <- sampleWeighted n $ map getweight minima
  let normalized = map (reweight $ Product (w / Exp (log $ fromIntegral n))) minima'
  return (normalized ++ rest)

resampleAlpha :: Int -> [AlphaSMC a] -> Random [AlphaSMC a]
resampleAlpha n ms = do
    -- do all the stochastic compuations so everything
    -- is either done or at a resample statement
    ms' <- parallel $ map (compute . unalpha) ms
    -- check if all particles are done...
    if all isDone ms'
      -- ... in which case return
      then return $ map AlphaSMC ms'
      -- otherwise we do resampling
      else do
        -- we need to prune the resample we dealt with and reset weights
        resampled <- resampleMinSubset ms'
        -- recursive call to continue inference
        resampleAlpha n $ map AlphaSMC resampled

resampleESS :: Int -> ([Adaptive a], Weight) -> Random ([Adaptive a], Weight)
resampleESS n (ms, w) = do
    -- do all the stochastic compuations so everything
    -- is either done or at a resample statement
    ms' <- parallel $ map (compute . unadapt) ms
    let w1 = sum (map weight ms')
    let w2 = sum (map ((^2) . weight) ms')
    let ess = w1 * w1 / w2
    let w' = w * (w1 / Exp (log $ fromIntegral n))

    -- check if all particles are done...
    if all isDone ms'
      -- ... in which case return
      then return (map Adaptive ms', w')
      -- otherwise we do resampling
      else do
        if ess < Exp (log (1 * fromIntegral n))
        then do
          resampled <- sampleWeighted n (map getweight ms')
          resampleESS n (map Adaptive resampled, w')
        else
          resampleESS n (map (Adaptive . ignoreResample) ms', w)

instance SMC AlphaSMC where
  resampleSMC n ms = do
    ms' <- resampleAlpha n ms
    return (map (\(AlphaSMC (Done w a)) -> (a, getProduct w)) ms')
  normalisingConstant n ms = do
    ms' <- resampleAlpha n ms
    let w = sum (map (\(AlphaSMC (Done w a)) -> getProduct w) ms')
    return $ w / Exp (log $ fromIntegral n)

instance SMC Labelled where
  resampleSMC n ms = do
    ms' <- resampleLab ms
    return (map (\(Labelled (Done w a)) -> (a, getProduct w)) ms')
  normalisingConstant n ms = do
    ms' <- resampleLab ms
    return $ getProduct (sum $ map (\(Labelled (Done w a)) -> w) ms') / Exp (log $ fromIntegral $ length ms')

instance SMC Adaptive where
  resampleSMC n ms = do
    (ms', w) <- resampleESS n (ms, Exp 0)
    return (map (\(Adaptive (Done w a)) -> (a, getProduct w)) ms')
  normalisingConstant n ms = snd <$> resampleESS n (ms, Exp 0)

instance SMC Stratified where
  resampleSMC n ms = do
    ms' <- resampleStrat ms
    return (map (\(Stratified (Done w a)) -> (a, getProduct w)) ms')
  normalisingConstant n ms = do
    ms' <- resampleStrat ms
    return $ getProduct (sum $ map (\(Stratified (Done w a)) -> w) ms') / Exp (log $ fromIntegral $ length ms')

ignoreResample :: ParticleT l (Product Weight) Random a -> ParticleT l (Product Weight) Random a
ignoreResample (Resample _ w a) = reweight w a
ignoreResample p@(Done _ _) = p
ignoreResample (Compute _) = error "unexpected compute"

instance Stochastic AliveFilter where
  -- keep trying until we get a nonzero weight
  compute p@(AliveFilter (Compute ma)) = do
    a <- compute (Compute ma)
    if weight a == Exp (log 0) -- note that log 0 = -Infinity
    then compute p
    else return $ AliveFilter a
  compute p = return p

instance Stochastic AlphaSMC where
  -- keep trying until we get a nonzero weight
  compute p@(AlphaSMC (Compute ma)) = do
    a <- compute (Compute ma)
    if weight a == Exp (log 0) -- note that log 0 = -Infinity
    then compute p
    else return $ AlphaSMC a
  compute p = return p

-- get weights, clear internal weights, remove resample
getweight :: ParticleT l (Product Weight) m a -> Weighted (ParticleT l (Product Weight) m a)
getweight (Done w a) = (Done mempty a, getProduct w)
getweight (Resample l w a) = (a, getProduct w)
getweight (Compute a) = error "unexpected compuation"

-- check if particle is done
isDone :: ParticleT l w m a -> Bool
isDone (Done _ _) = True
isDone _ = False

smcmodel :: Int -> ([SMCModel a], Weight) -> Random ([SMCModel a], Weight)
smcmodel n (ms, w) = do
  -- do all the stochastic compuations so everything
    -- is either done or at a resample statement
    ms' <- parallel $ map compute ms
    let w' = w * sum (map weight ms') / Exp (log $ fromIntegral n)
    -- check if all particles are done...
    if all isDone ms'
      -- ... in which case return
      then return (ms', w')
      -- otherwise we do resampling
      else do
        -- we need to prune the resample we dealt with and reset weights
        resampled <- sampleWeighted n (map getweight ms')
        -- recursive call to continue inference
        smcmodel n (resampled, w')

instance SMC SMCModel where
  resampleSMC n ms = do
    (ms', w) <- smcmodel n (ms, Exp 0)
    return (map (\(Done w a) -> (a, getProduct w)) ms')
  normalisingConstant n ms = snd <$> smcmodel n (ms, Exp 0)

inferSMC :: SMC m => Int -> m a -> Random [Weighted a]
inferSMC n m = resampleSMC n (replicate n m)

inferConstant :: SMC m => Int -> m a -> Random Weight
inferConstant n m = normalisingConstant n (replicate n m)

max_stars_repo_name: kai-pischke/ppl

max_issues_repo_name: kai-pischke/ppl

max_forks_repo_name: kai-pischke/ppl


text:
-- The magic of Lie Algebra

import Numeric.LinearAlgebra

disp = putStrLn . dispf 5

rot1 :: Double -> Matrix Double
rot1 a = (3><3)
 [ 1, 0, 0
 , 0, c, s
 , 0,-s, c ]
    where c = cos a
          s = sin a

g1,g2,g3 :: Matrix Double

g1 = (3><3) [0, 0,0
            ,0, 0,1
            ,0,-1,0]

rot2 :: Double -> Matrix Double
rot2 a = (3><3)
 [ c, 0, s
 , 0, 1, 0
 ,-s, 0, c ]
    where c = cos a
          s = sin a

g2 = (3><3) [ 0,0,1
            , 0,0,0
            ,-1,0,0]

rot3 :: Double -> Matrix Double
rot3 a = (3><3)
 [ c, s, 0
 ,-s, c, 0
 , 0, 0, 1 ]
    where c = cos a
          s = sin a

g3 = (3><3) [ 0,1,0
            ,-1,0,0
            , 0,0,0]

deg=pi/180

-- commutator
infix 8 &
a & b = a <> b - b <> a

infixl 6 |+|
a |+| b = a + b + a&b /2  + (a-b)&(a & b) /12

main = do
   let a =  45*deg
       b =  50*deg
       c = -30*deg
       exact = rot3 a <> rot1 b <> rot2 c
       lie = scalar a * g3 |+| scalar b * g1 |+| scalar c * g2
   putStrLn "position in the tangent space:"
   disp lie
   putStrLn "exponential map back to the group (2 terms):"
   disp (expm lie)
   putStrLn "exact position:"
   disp exact

max_stars_repo_name: curiousleo/liquidhaskell

max_issues_repo_name: curiousleo/liquidhaskell

max_forks_repo_name: curiousleo/liquidhaskell


text:
-- |
-- Module      : Main
-- Description :
-- Copyright   : (c) Jonatan H Sundqvist, 2015
-- License     : MIT
-- Maintainer  : Jonatan H Sundqvist
-- Stability   : experimental|stable
-- Portability : POSIX (not sure)
--

-- Created September 20 2015

-- TODO | -
--        -

-- SPEC | -
--        -

--------------------------------------------------------------------------------------------------------------------------------------------
-- GHC Pragmas
--------------------------------------------------------------------------------------------------------------------------------------------

--------------------------------------------------------------------------------------------------------------------------------------------
-- API
--------------------------------------------------------------------------------------------------------------------------------------------
module Main where

--------------------------------------------------------------------------------------------------------------------------------------------
-- We'll need these
--------------------------------------------------------------------------------------------------------------------------------------------
import Data.IORef
import Data.Complex
import qualified Data.Set as S
import qualified Data.Map as M

import Graphics.UI.Gtk
import qualified Graphics.Rendering.Cairo as Cairo

import Southpaw.Math.Constants
import qualified Southpaw.Picasso.Palette as Palette

import Occlusion.Types
import Occlusion.Vector
import qualified Occlusion.Window     as Window
import qualified Occlusion.Events     as Events
import qualified Occlusion.Behaviours as Behaviours

--------------------------------------------------------------------------------------------------------------------------------------------
-- Entry point
--------------------------------------------------------------------------------------------------------------------------------------------
-- |
main :: IO ()
main = do
  (window, canvas) <- Window.create size
  assets'          <- loadAssets
  stateref         <- newIORef $ AppState { _animation = AnimationData { _fps=30.0, _frame=0 },
                                            _input     = InputData { _keyboard=S.empty, _mouse=0:+0, _click=Nothing },
                                            _scene     = Scene { _player=Character { _position=0:+0, _velocity=0:+0, _health=100, _name="Democritus", _colour=Palette.peru,  _behaviour=Behaviours.player },
                                                                 _npcs=[Character { _position=0:+0, _velocity=20:+20, _health=58, _name="Agamemnon", _colour=Palette.orchid, _behaviour=Behaviours.npc patrol }],
                                                                 _obstacles=theobstacles },
                                            _gui       = GUI { _window=window, _canvas=canvas },
                                            _assets    = assets' }
  Events.attach window canvas stateref
  mainGUI
  where
    size = (720:+480) * 1.2
    patrol = map (+(200:+150)) [20:+20, 300:+32, 200:+10, 20:+20]

-- |
loadAssets :: IO Assets
loadAssets = do
  tree <- Cairo.imageSurfaceCreateFromPNG "assets/images/TreeAtDuskSmall.png"
  return $ Assets { _images=M.fromList [("tree", tree)] }

-- |
theobstacles :: [Polygon Double]
theobstacles = [[origin + let θ = α + n*2*π/ns in mkPolar mag θ | n <- [0..(ns-1)]] | (α, ns, mag, origin) <- [(0, 5, 78, 522:+150), (0, 8, 92, 192:+342)] ]

max_stars_repo_name: SwiftsNamesake/Occlusion

max_issues_repo_name: SwiftsNamesake/Occlusion

max_forks_repo_name: SwiftsNamesake/Occlusion


text:
{-# LANGUAGE FlexibleContexts #-}
module PlotGreensFunctionSTS where

import           Control.Monad.Parallel    as MP
import           Data.Array.Repa           as R
import           Data.Binary               (decodeFile)
import           Data.Complex
import           Data.List                 as L
import           Data.Vector.Unboxed       as VU
import           FokkerPlanck.DomainChange
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.IO
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types
import           Utils.Array

main = do
  args@(numPointStr:spatialFreqStr:numOrientationStr:thetaSigmaStr:numScaleStr:scaleSigmaStr:maxScaleStr:taoStr:initStr:numTrailStr:maxTrailStr:thetaFreqsStr:scaleFreqsStr:histFileName:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      spatialFreq = read spatialFreqStr :: Double
      spatialFreqs = [-spatialFreq .. spatialFreq]
      numOrientation = read numOrientationStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      numScale = read numScaleStr :: Int
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      init@(initX, initY, initTheta, initScale) =
        read initStr :: (Double, Double, Double, Double)
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      thetaFreq = read thetaFreqsStr :: Double
      scaleFreq = read scaleFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/PlotGreensFunctionSTS"
      histPath = folderPath </> histFileName
      initArr =
        traverse4
          (fromListUnboxed (Z :. (L.length spatialFreqs)) spatialFreqs)
          (fromListUnboxed (Z :. (L.length spatialFreqs)) spatialFreqs)
          (fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs)
          (fromListUnboxed (Z :. (L.length scaleFreqs)) scaleFreqs)
          (\(Z :. a) (Z :. b) (Z :. c) (Z :. d) ->
             (Z :. a :. b :. c :. d :. c :. d)) $ \fx fy ft fs (Z :. x :. y :. _ :. _ :. t :. s) ->
          exp
            (0 :+
             (-1) *
             (ft (Z :. t) * (initTheta / 180 * pi) +
              2 * pi *
              ((fx (Z :. x) * initX + fy (Z :. y) * initY) /
               (fromIntegral numPoint) +
               fs (Z :. s) * initScale / log maxScale)))
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histPath
  arr <-
    if flag
      then fmap
             (fromListUnboxed
                (Z :. (L.length spatialFreqs) :. (L.length spatialFreqs) :.
                 (L.length thetaFreqs) :.
                 (L.length scaleFreqs) :.
                 (L.length thetaFreqs) :.
                 (L.length scaleFreqs))) .
           decodeFile $
           histPath
      else runMonteCarloSTS
             numThread
             numTrail
             maxTrail
             numPoint
             numPoint
             thetaSigma
             scaleSigma
             maxScale
             tao
             spatialFreqs
             spatialFreqs
             thetaFreqs
             scaleFreqs
             histPath $
           VU.replicate
             (L.length scaleFreqs * L.length thetaFreqs * L.length scaleFreqs *
              L.length thetaFreqs *
              L.length spatialFreqs *
              L.length spatialFreqs)
             0
  transformedArr <- sumP . sumS $ R.zipWith (*) arr initArr
  -- arrR2S1RP <-
  --   stsTor2s1rp
  --     numPoint
  --     (fromIntegral numPoint)
  --     spatialFreqs
  --     numPoint
  --     (fromIntegral numPoint)
  --     spatialFreqs
  --     numOrientation
  --     thetaFreqs
  --     numScale
  --     scaleFreqs
  --     maxScale
  --     transformedArr
  -- plotImageRepaComplex (folderPath </> "Greens.png") .
  --   ImageRepa 8 .
  --   computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.sumS . sumS $
  --   arrR2S1RP
  arrMag <-
    stsTor2'
      numPoint
      (fromIntegral numPoint)
      spatialFreqs
      numPoint
      (fromIntegral numPoint)
      spatialFreqs
      transformedArr
  plotImageRepa (folderPath </> "GreensMag.png") .
    ImageRepa 8 .
    computeS . R.extend (Z :. (1 :: Int) :. All :. All) . R.map sqrt $
    arrMag

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field