text:
module PropertiesTests
 ( properties
 ) where

import GradientDescent
import Test.Tasty
import qualified Test.Tasty.QuickCheck as QC
import qualified Test.Tasty.SmallCheck as SC
import Numeric.LinearAlgebra

properties :: TestTree
properties = testGroup "Properties" [scProps, qcProps]

scProps =
  testGroup
    "(checked by SmallCheck)"
    [checkTranspositionSC, checkMatrixSigmoidSC, checkMatrixSigmoidPrimeSC]

checkTranspositionSC = SC.testProperty "matrix == tr . tr matrix" $
                           \list -> ((length list >< 1) (list :: [Double]) :: Matrix Double)
                              == (tr $ tr $ (length list><1) (list :: [Double]) :: Matrix Double)

checkMatrixSigmoidSC = SC.testProperty "check matrixSigmoid" $
                         \list ->
                            (length (list :: [Double])) >=
                              1 SC.==> ((1 >< length list) (sigmoid <$> list) :: Matrix Double)
                                == matrixSigmoid ((1 >< length list) list :: Matrix Double)

checkMatrixSigmoidPrimeSC = SC.testProperty "check matrixSigmoidPrime" $
                         \list ->
                            (length (list :: [Double])) >=
                              1 SC.==> ((1 >< length list) (sigmoidPrime <$> list) :: Matrix Double)
                                == matrixSigmoidPrime ((1 >< length list) list :: Matrix Double)

qcProps =
  testGroup
    "(checked by QuickCheck)"
    [checkTranspositionQC, checkMatrixSigmoidQC, checkMatrixSigmoidPrimeQC]

checkTranspositionQC = QC.testProperty "matrix == tr . tr matrix" $
                           \list -> ((length list >< 1) (list :: [Double]) :: Matrix Double)
                              == (tr $ tr $ (length list><1) (list :: [Double]) :: Matrix Double)

checkMatrixSigmoidQC = QC.testProperty "check matrixSigmoid" $
                         \list ->
                            (length (list :: [Double])) >=
                              1 QC.==> ((1 >< length list) (sigmoid <$> list) :: Matrix Double)
                                == matrixSigmoid ((1 >< length list) list :: Matrix Double)

checkMatrixSigmoidPrimeQC = QC.testProperty "check matrixSigmoidPrime" $
                         \list ->
                            (length (list :: [Double])) >=
                              1 QC.==> ((1 >< length list) (sigmoidPrime <$> list) :: Matrix Double)
                                == matrixSigmoidPrime ((1 >< length list) list :: Matrix Double)

max_stars_repo_name: Malenczuk/HuskNet

max_issues_repo_name: Malenczuk/HuskNet

max_forks_repo_name: Malenczuk/HuskNet


text:
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util (ones,zeros)
import HVX

n = 4
a = EConst $ (n><n) ([0, 2, 0, 0,
                      -3, 0, 0, 0,
                       0, 0, 1, 0,
                       0, 0, 0, 4] :: [Double])
b = EConst $ (n><n) ([4, 8, -1, 0,
                      0, -2, 1, 0,
                      9, -3, 0, 3,
                      0, 0,  2, 0] :: [Double])
c = EConst $ (n><1) [1..]

x = EVar "x"
y = EVar "y"
constNeg3 = EConst $ scale (-3.0) $ ones n 1
constTwos = EConst $ scale 2.0 $ ones n 1

ans = subgradMinimize
  (norm 2 (a *~ x) +~ norm 2 (b *~ y) +~ norm 2 (c +~ x +~ neg y))
  [y <=~ constTwos, x <=~ constNeg3]
  (decNonSumStep 10.0) 5000
  [("x",zeros n 1), ("y", zeros n 1)]

main :: IO ()
main = do
  putStrLn $ show ans

ext: hs

max_stars_repo_name: wellposed/hvx

max_issues_repo_name: wellposed/hvx

max_forks_repo_name: wellposed/hvx


text:
{-|
Module: MachineLearning.Classification.OneVsAll
Description: One-vs-All Classification.
Copyright: (c) Alexander Ignatyev, 2016-2017
License: BSD-3
Stability: experimental
Portability: POSIX

One-vs-All Classification.
-}

module MachineLearning.Classification.OneVsAll
(
    Opt.MinimizeMethod(..)
  , module Log
  , module Model
  , predict
  , learn
  , MLC.calcAccuracy
  , Regularization(..)
)

where

import MachineLearning.Types (R, Vector, Matrix)
import MachineLearning.Regularization (Regularization(..))
import qualified MachineLearning.Optimization as Opt
import qualified MachineLearning.LogisticModel as Log
import qualified MachineLearning.Model as Model
import qualified MachineLearning.Classification.Internal as MLC
import qualified Data.Vector.Storable as V
import qualified Numeric.LinearAlgebra as LA

-- | One-vs-All Classification prediction function.
-- Takes a matrix of features X and a list of vectors theta,
-- returns predicted class number assuming that class numbers start at 0.
predict :: Matrix -> [Vector] -> Vector
predict x thetas = predictions'
  where predict = Model.hypothesis Log.Logistic x
        predictions = LA.toRows . LA.fromColumns $ map predict thetas
        predictions' = LA.vector $ map (fromIntegral . LA.maxIndex) predictions

-- | Learns One-vs-All Classification
learn :: Opt.MinimizeMethod -- ^ (e.g. BFGS2 0.1 0.1)
         -> R                  -- ^ epsilon, desired precision of the solution;
         -> Int                -- ^ maximum number of iterations allowed;
         -> Regularization     -- ^ regularization parameter lambda;
         -> Int                -- ^ number of labels
         -> Matrix             -- ^ matrix X;
         -> Vector             -- ^ vector y
         -> [Vector]             -- ^ initial theta list;
         -> ([Vector], [Matrix])  -- ^ solution vector and optimization path.
learn mm eps numIters lambda nLabels x y initialThetaList =
  let ys = MLC.processOutputOneVsAll nLabels y
      minimize = Opt.minimize mm Log.Logistic eps numIters lambda x
  in unzip $ zipWith minimize ys initialThetaList

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module AutoEncoder (
  preTrains,
  preTrain
  ) where

import Numeric.LinearAlgebra
import Data.List
import Common
import Forward
import BackProp
import ActivFunc
import Other

-- TODO: refactor
preTrains :: R -> Int -> (Matrix R -> Matrix R, Matrix R -> Matrix R) -> Matrix R -> [Matrix R] -> [Matrix R]
preTrains rate iter fDf x ws = snd (mapAccumL (preTrain rate iter fDf) x ms) `mappend` [m]
  where (ms, m) = (init ws, last ws)

preTrain :: R -> Int -> (Matrix R -> Matrix R, Matrix R -> Matrix R) -> Matrix R -> Matrix R -> (Matrix R, Matrix R)
preTrain rate iter fDf@(f, df) x w = (nx, nw)
  where
    nx = forward f nw x
    [nw, _] = last . take iter $ iterate (backProp f rate fDf (x, x)) [w, mirrorWeight w]

mirrorWeight :: Matrix R -> Matrix R
mirrorWeight w = w' ||| konst 0 (rows w', 1) -- initial bias: 0
  where w' = tr $ weightWithoutBias w

max_stars_repo_name: pupuu/deep-neuralnet

max_issues_repo_name: pupuu/deep-neuralnet

max_forks_repo_name: pupuu/deep-neuralnet


text:
module Section_3_4_3 where

import Numeric.LinearAlgebra

type Network = [(Matrix Double, Vector Double, Double -> Double)]

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

initNetwork :: Network
initNetwork =
  [ ( (2><3) [0.1, 0.3, 0.5, 0.2, 0.4, 0.6]
    , vector [0.1, 0.2, 0.3]
    , sigmoid
    )
  , ( (3><2) [0.1, 0.4, 0.2, 0.5, 0.3, 0.6]
    , vector [0.1, 0.2]
    , sigmoid
    )
  , ( (2><2) [0.1, 0.3, 0.2, 0.4]
    , vector [0.1, 0.2]
    , id
    )
  ]

forward :: Network -> Vector Double -> Vector Double
forward nw x = foldl (\x (w, b, h) -> cmap h (x <# w + b)) x nw


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

import Criterion.Config
import Criterion.Main
import Statistics.Distribution.Normal

import Statistics.Distribution.Normal
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU
import qualified Data.Vector.Generic as G

import Data.Array.Unboxed

import HLearn.Algebra
import HLearn.Models.Distributions
import HLearn.Models.Distributions.Moments
import HLearn.Models.Distributions.MultiNormal
-- import qualified HLearn.Models.Distributions.GaussianOld as GO
-- import qualified HLearn.Models.Distributions.GaussianOld2 as GO2

import qualified Control.ConstraintKinds as CK

-- comparing different types

{-size = 10^6
main = defaultMain
    [ bench "Moments Double" $ nf (train :: [Double] -> Moments Double) [(0::Double)..size]
    ]-}

myConfig = defaultConfig
    { cfgPerformGC = ljust True
    , cfgSamples = ljust 5
    }

size = 10^6
main = defaultMainWith myConfig (return ())
-- main = defaultMain
    [ bench "HLearn-MultiNormalVec 1" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 1 Double) $ map VU.fromList [ [l] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 2" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 2 Double) $ map VU.fromList [ [l..l+1] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 3" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 3 Double) $ map VU.fromList [ [l..l+2] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 4" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 4 Double) $ map VU.fromList [ [l..l+3] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 5" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 5 Double) $ map VU.fromList [ [l..l+4] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 10" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 10 Double) $ map VU.fromList [ [l..l+9] | l <- [(0::Double)..fromIntegral size]]
    , bench "HLearn-MultiNormalVec 100" $ whnf (train :: [VU.Vector Double] -> MultiNormalVec 100 Double) $ map VU.fromList [ [l..l+99] | l <- [(0::Double)..fromIntegral size]]
    ]
-- [ bench "HLearn-MultiNormal-Vector" $ whnf
--         ((train :: [VU.Vector Double] -> MultiNormalVec 1 Double))
--         $ map VU.singleton [(0::Double)..fromIntegral size]
--     , bench "HLearn-MultiNormal" $ whnf
--         ((train :: [UArray Int Double] -> MultiNormal Double 1))
--         $ map (\x -> listArray (0,0) [x]) [(0::Double)..fromIntegral size]
--     ]

{-main = defaultMain
    [ bench "HLearn-GaussianCK" $ whnf
        (trainCK :: V.Vector Double -> Gaussian Double)
--         [(0::Double)..fromIntegral size]
        (V.enumFromN (0::Double) size)
    , bench "HLearn-Gaussian" $ whnf
        (train :: V.Vector Double -> Gaussian Double)
--         [(0::Double)..fromIntegral size]
        (V.enumFromN (0::Double) size)
    , bench "HLearn-Moments" $ whnf
        ((train :: V.Vector Double -> Moments Double))
        (V.enumFromN (0::Double) size)
    , bench "HLearn-MomentsSingleton" $ whnf
        ((train :: [VU.Vector Double] -> NormSingle Double))
        $ map VU.singleton [(0::Double)..fromIntegral size]
    , bench "HLearn-MultiNormal-Vector" $ whnf
        ((train :: [VU.Vector Double] -> MultiNormal2 1 Double))
        $ map VU.singleton [(0::Double)..fromIntegral size]
    , bench "HLearn-MultiNormal" $ whnf
        ((train :: [UArray Int Double] -> MultiNormal Double 1))
        $ map (\x -> listArray (0,0) [x]) [(0::Double)..fromIntegral size]
    ]

newtype NormSingle a = NormSingle (Moments a)
    deriving (Read,Show,Eq,Ord,Semigroup,Monoid)

instance Model (NoParams (NormSingle a)) (NormSingle a) where
    getparams _ = NoParams

instance DefaultModel (NoParams (NormSingle a)) (NormSingle a) where
    defparams = NoParams

instance (Num a, VU.Unbox a) => HomTrainer (NoParams (NormSingle a)) (VU.Vector a) (NormSingle a) where
    train1dp' params dp = NormSingle $ train1dp (dp VU.! 0)-}

-- size = 10^6
-- main = defaultMain
--     [ bench "Moments 2 Double" $ nf (train :: V.Vector Double -> Moments Double) (V.enumFromN (0::Double) size)
--     , bench "Moments 2 Float" $ nf (train :: V.Vector Float -> Moments Float) (V.enumFromN (0::Float) size)
--     , bench "Moments 2 Rational" $ nf (train :: V.Vector Rational -> Moments Rational) (V.enumFromN (0::Rational) size)
--     ]

{-main = defaultMain
    [ bench "HLearn-Gaussian" $ nf ((train :: VU.Vector Double -> Gaussian Double)) (VU.enumFromN (0::Double) size)
    , bench "HLearn-Gaussian-Parallel" $ whnf (parallel $ (train :: VU.Vector Double -> Gaussian Double)) (VU.enumFromN (0::Double) size)
--     , bench "HLearn-Gaussian-List" $ nf (train :: [Double] -> Gaussian Double) [0..fromIntegral size]
    , bench "statistics-Gaussian" $ whnf (normalFromSample . VU.enumFromN 0) (size)
    ]-}

    {-bench "batch train [] 1e6" $ nf ((batch train) GaussianParams) [0..1e6::Double]
--     , bench "batch train V 1e6"  $ nf ((batch train) GaussianParams) (V.enumFromN (0::Double) (10^6))
--     , bench "batch train VU 1e6" $ nf ((batch train) GaussianParams) (VU.enumFromN (0::Double) [0..1e6::Double])
--     , bench "parallel2 batch train [] 1e6" $ nf ((parallel $ batch train) GaussianParams) (V.enumFromN (0::Double) (10^7))
    ,-}
--     bench "Parallel GaussianOld2" $ nf (parallel $ batch (train GO2.GaussianParams)) (VU.enumFromN (0::Double) (10^8))
--     , bench "GaussianOld2" $ nf (batch (train GO2.GaussianParams)) (VU.enumFromN (0::Double) (10^8))
--     bench "Parallel GaussianOld" $ nf (parallel $ batch (train GO.GaussianParams)) (VU.enumFromN (0::Double) (10^8))
--     , bench "GaussianOld" $ nf (batch (train GO.GaussianParams)) (VU.enumFromN (0::Double) (10^8))
--      bench "Parallel Gaussian2" $ nf (parallel $ batch (train GaussianParams)) (VU.enumFromN (0::Double) (10^8))
--     , bench "Gaussian2" $ nf (batch (train GaussianParams)) (VU.enumFromN (0::Double) (10^8))

    {-bench "Parallel Gaussian" $ nf (parallel $ batch (trainSG)) (VU.enumFromN (0::Double) (10^8))
    , bench "Gaussian" $ nf (batch (trainSG)) (VU.enumFromN (0::Double) (10^8))
    ,-} -- bench "Categorical" $ nf (train :: [String] -> Categorical String Double) (concat $ replicate (10^4) ["a","b"])
--     , bench "normalFromSample - 1e5" $ whnf (normalFromSample . VU.enumFromN 1) (10^8)

-- instance NFData NormalDistribution where
--     rnf (NormalDistribution a b c d) = deepseq a $ deepseq b $ deepseq c $ rnf d

-- parallel3 strat train = \modelparams datapoint ->
--     foldl1 (<>) $ map (CK.foldl1 (<>) . CK.fmap (train modelparams . CK.pure)) (CK.partition 2 datapoint)

-- test (ds1,ds2) = m1 <> m2
--     where
--         [m1,m2] = parMap rdeepseq ((batch train) GaussianParams) [ds1,ds2]
-- --        [m1,m2] = parMap rdeepseq ((batch train) GaussianParams) [ds1,ds2]
--
-- -- test (ds1,ds2) = runEval $ do
-- --     let [m1,m2] = {-parMap rdeepseq-} map (batch train GaussianParams) [ds1,ds2]
-- -- --     m1 <- rparWith rseq $ batch train GaussianParams ds1
-- -- --     m2 <- rparWith rseq $ batch train GaussianParams ds2
-- -- {-    let m1 = batch train GaussianParams ds1
-- --     let m2 = batch train GaussianParams ds2-}
-- --     return $ m1 <> m2
--
-- dsL = replicate 4 ds
-- -- dsL = [ds1,ds2]
-- ds = [0..1000000::Double]
-- ds1 = V.enumFromN (0::Double) (10^7)
-- ds2 = V.enumFromN (0.5::Double) (10^7)

max_stars_repo_name: silky/HLearn


text:
{-# LANGUAGE DeriveGeneric    #-}
{-# LANGUAGE TypeApplications #-}
module Noether.Algebra.Actions.Acts where

import           Data.Complex

import           Noether.Lemmata.Prelude
import           Noether.Lemmata.TypeFu

import           Noether.Algebra.Single
import           Noether.Algebra.Tags

import           GHC.Generics

data ActsE
  = Acts_Magma MagmaE
  | ActsNamed Symbol ActsE
  | ActsTagged Type ActsE
  deriving (Generic)

class ActsK (lr :: Side) (op :: k) a b (s :: ActsE) where
  actK :: Proxy op -> Proxy s -> Proxy lr -> a -> b -> b

type family ActsS (lr :: Side) (op :: k) (a :: Type) (b :: Type) = (r :: ActsE)

instance MagmaK op a zm => ActsK lr op a a (Acts_Magma zm) where
  actK opP _ _ = binaryOpK opP (Proxy @zm)

max_stars_repo_name: evertedsphere/noether

max_issues_repo_name: evertedsphere/noether

max_forks_repo_name: evertedsphere/noether


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE DeriveGeneric       #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE LambdaCase          #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving  #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

import Control.Monad.Random
import Data.Binary                  as B
import Data.Singletons
import Data.Singletons.Prelude
import Data.Singletons.TypeLits
import GHC.Generics                 (Generic)
import Numeric.LinearAlgebra.Static

-- * Weights
--

data Weights i o = W { wBiases :: !(R o)
                     , wNodes  :: !(L o i)
                     }
  deriving (Show, Generic)

instance (KnownNat i, KnownNat o) => Binary (Weights i o)

randomWeights :: (MonadRandom m, KnownNat i, KnownNat o)
              => m (Weights i o)
randomWeights = do
    s1 :: Int <- getRandom
    s2 :: Int <- getRandom
    let wB = randomVector  s1 Uniform * 2 - 1
        wN = uniformSample s2 (-1) 1
    return $ W wB wN

runLayer :: (KnownNat i, KnownNat o)
         => Weights i o
         -> R i
         -> R o
runLayer (W wB wN) v = wB + wN #> v

-- * Network
--

data Network :: Nat -> [Nat] -> Nat -> * where
    O     :: !(Weights i o)
          -> Network i '[] o
    (:&~) :: KnownNat h
          => !(Weights i h)
          -> !(Network h hs o)
          -> Network i (h ': hs) o
infixr 5 :&~

deriving instance (KnownNat i, KnownNat o) => Show (Network i hs o)

logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))

runNet :: (KnownNat i, KnownNat o)
       => Network i hs o
       -> R i
       -> R o
runNet (O w)      !v = logistic (runLayer w v)
runNet (w :&~ n') !v = let v' = logistic (runLayer w v)
                       in  runNet n' v'

hiddenStruct :: Network i hs o -> [Integer]
hiddenStruct = \case
    O _    -> []
    _ :&~ (n' :: Network h hs' o)
           -> natVal (Proxy @h)
            : hiddenStruct n'

randomNet' :: forall m i hs o. (MonadRandom m, KnownNat i, KnownNat o)
           => Sing hs -> m (Network i hs o)
randomNet' = \case
    SNil            ->     O <$> randomWeights
    SNat `SCons` ss -> (:&~) <$> randomWeights <*> randomNet' ss

randomNet :: forall m i hs o. (MonadRandom m, KnownNat i, SingI hs, KnownNat o)
          => m (Network i hs o)
randomNet = randomNet' sing

putNet :: (KnownNat i, KnownNat o)
       => Network i hs o
       -> Put
putNet = \case
    O w     -> put w
    w :&~ n -> put w *> putNet n

getNet :: forall i hs o. (KnownNat i, KnownNat o)
       => Sing hs
       -> Get (Network i hs o)
getNet = \case
    SNil            ->     O <$> get
    SNat `SCons` ss -> (:&~) <$> get <*> getNet ss

instance (KnownNat i, SingI hs, KnownNat o) => Binary (Network i hs o) where
    put = putNet
    get = getNet sing

-- * OpaqueNet
--

data OpaqueNet :: Nat -> Nat -> * where
    ONet :: Network i hs o -> OpaqueNet i o

runOpaqueNet :: (KnownNat i, KnownNat o)
             => OpaqueNet i o
             -> R i
             -> R o
runOpaqueNet (ONet n) x = runNet n x

numHiddens :: OpaqueNet i o -> Int
numHiddens (ONet n) = go n
  where
    go :: Network i hs o -> Int
    go = \case
        O _      -> 0
        _ :&~ n' -> 1 + go n'

randomONet :: (MonadRandom m, KnownNat i, KnownNat o)
           => [Integer]
           -> m (OpaqueNet i o)
randomONet hs = case toSing hs of
                  SomeSing ss -> ONet <$> randomNet' ss

putONet :: (KnownNat i, KnownNat o)
        => OpaqueNet i o
        -> Put
putONet (ONet net) = do
    put (hiddenStruct net)
    putNet net

getONet :: (KnownNat i, KnownNat o)
        => Get (OpaqueNet i o)
getONet = do
    hs <- get
    withSomeSing hs $ \ss ->
      ONet <$> getNet ss

instance (KnownNat i, KnownNat o) => Binary (OpaqueNet i o) where
    put = putONet
    get = getONet

-- * OpaqueNet'
--

type OpaqueNet' i o r = (forall hs. Network i hs o -> r) -> r

oNet' :: Network i hs o -> OpaqueNet' i o r
oNet' n = \f -> f n

runOpaqueNet' :: (KnownNat i, KnownNat o)
              => OpaqueNet' i o (R o)
              -> R i
              -> R o
runOpaqueNet' oN x = oN (\n -> runNet n x)
--            :: ((forall hs. Network i hs o -> R o) -> R o)
--            -> R i
--            -> R o

numHiddens' :: OpaqueNet' i o Int -> Int
numHiddens' oN = oN go
  where
    go :: Network i hs o -> Int
    go = \case
        O _      -> 0
        _ :&~ n' -> 1 + go n'
--          :: ((forall hs. Network i hs o -> Int) -> Int)
--          -> Int

withRandomONet' :: (MonadRandom m, KnownNat i, KnownNat o)
                => [Integer]
                -> (forall hs. Network i hs o -> m r)
                -> m r
--         aka, => [Integer]
--              -> OpaqueNet' i o (m r)
withRandomONet' hs f = withSomeSing hs $ \ss -> do
                         net <- randomNet' ss
                         f net

putONet' :: (KnownNat i, KnownNat o)
         => OpaqueNet' i o Put
         -> Put
putONet' oN = oN $ \net -> do
                      put (hiddenStruct net)
                      putNet net

getONet' :: (KnownNat i, KnownNat o)
         => (forall hs. Network i hs o -> Get r)
         -> Get r
--  aka, => OpaqueNet' i o (Get r)
getONet' f = do
    hs <- get :: Get [Integer]
    withSomeSing hs $ \ss -> do
      n <- getNet ss
      f n

main :: IO ()
main = do
    putStrLn "What hidden layer structure do you want?"
    hs <- readLn
    n  <- randomONet hs
    case n of
      ONet (net :: Network 10 hs 3) -> do
        print net
        -- blah blah stuff with our dynamically generated net

main' :: IO ()
main' = do
    putStrLn "What hidden layer structure do you want?"
    hs <- readLn
    withRandomONet' hs $ \(net :: Network 10 hs 3) -> do
      print net
      -- blah blah stuff with our dynamically generated net

-- * SomeNet
--

data SomeNet where
    SNet :: (KnownNat i, KnownNat o)
         => Network i hs o
         -> SomeNet

withONet :: SomeNet
         -> (forall i o. (KnownNat i, KnownNat o) => OpaqueNet i o -> r)
         -> r
withONet (SNet n) f = f (ONet n)

randomSNet :: forall m. MonadRandom m
           => Integer
           -> [Integer]
           -> Integer
           -> m SomeNet
randomSNet i hs o =
    withSomeSing i  $ \(SNat :: Sing (i  :: Nat  )) ->
    withSomeSing hs $ \(ss   :: Sing (hs :: [Nat])) ->
    withSomeSing o  $ \(SNat :: Sing (o  :: Nat  )) ->
      SNet <$> (randomNet' ss :: m (Network i hs o))

withRandomSNet' :: forall m r. MonadRandom m
                => Integer
                -> [Integer]
                -> Integer
                -> (forall i hs o. (KnownNat i, KnownNat o) => Network i hs o -> m r)
                -> m r
withRandomSNet' i hs o f =
    withSomeSing i  $ \(SNat :: Sing (i  :: Nat  )) ->
    withSomeSing hs $ \(ss   :: Sing (hs :: [Nat])) ->
    withSomeSing o  $ \(SNat :: Sing (o  :: Nat  )) -> do
      n <- randomNet' ss :: m (Network i hs o)
      f n

instance Binary SomeNet where
    put (SNet (net :: Network i hs o)) = do
        put $ natVal (Proxy @i)
        put $ hiddenStruct net
        put $ natVal (Proxy @o)
        putNet net
    get = do
        i  <- get :: Get Integer
        hs <- get :: Get [Integer]
        o  <- get :: Get Integer
        withSomeSing i    $ \(SNat :: Sing (i  :: Nat  )) ->
          withSomeSing hs $ \(ss   :: Sing (hs :: [Nat])) ->
          withSomeSing o  $ \(SNat :: Sing (o  :: Nat  )) -> do
            n <- getNet ss :: Get (Network i hs o)
            return (SNet n)

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
{-# LANGUAGE BangPatterns #-}
module STCPinwheel where

import           Control.Monad                  as M
import           Control.Monad.Parallel         as MP
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.IArray              as IA
import           Data.Array.Repa                as R
import           Data.Binary
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           Filter.Utils
import           FokkerPlanck
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2D
import           FourierPinwheel
import           Image.IO
import           Pinwheel.FourierSeries2D
import           STC                            hiding (convolve)
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time

main = do
  args@(deviceIDsStr:numPointsStr:deltaStr:thresholdStr:numPointsReconStr:deltaReconStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:tauStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:stdR2Str:stdStr:numBatchR2Str:numBatchR2FreqsStr:numBatchOriStr:batchSizeStr:sStr:radiusStr:deltaTStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      delta = read deltaStr :: Double
      threshold = read thresholdStr :: Double
      numPointsRecon = read numPointsReconStr :: Int
      deltaRecon = read deltaReconStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tau = read tauStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      initSource = [L.head initPoints]
      initSink = [L.last initPoints]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCPinwheel"
      stdR2 = read stdR2Str :: Double
      std = read stdStr :: Double
      numBatchR2 = read numBatchR2Str :: Int
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      numBatchOri = read numBatchOriStr :: Int
      batchSize = read batchSizeStr :: Int
      s = read sStr :: Double
      radius = read radiusStr :: Double
      deltaT = read deltaTStr :: Double
      periodEnvelope = periodR2 ^ 2 * 2
      maxScale = sqrt periodEnvelope
  removePathForcibly folderPath
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath

  hist <-
    if flag
      then do printCurrentTime "read from files..."
              decodeFile histFilePath
      else do printCurrentTime "Start computing coefficients..."
              initialise []
              devs <- M.mapM device deviceIDs
              ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
              ptxs <- M.mapM createTargetFromContext ctxs
              runMonteCarloFourierCoefficientsGPU
                deviceIDs
                numThread
                1000000
                500000
                thetaSigma
                scaleSigma
                0
                s
                tau
                deltaT
                phiFreqs
                rhoFreqs
                thetaFreqs
                scaleFreqs
                maxScale
                histFilePath
      -- else sampleCartesian
      --        histFilePath
      --        folderPath
      --        ptxs
      --        numPoints
      --        periodEnvelope
      --        delta
      --        numOrientation
      --        initScale
      --        thetaSigma
      --        tau
      --        threshold
      --        s
      --        phiFreq
      --        rhoFreq
      --        thetaFreq
      --        scaleFreq
  printCurrentTime "Done"
  printCurrentTime "Start Convloution.."
  plan <-
    makePlanDiscrete
      folderPath
      emptyPlan
      numPointsRecon
      numPointsRecon
      numR2Freq
      (2 * thetaFreq + 1)
      (2 * scaleFreq + 1)
      (2 * phiFreq + 1)
      (2 * rhoFreq + 1)
  let coefficients = getNormalizedHistogramArr $ hist
  -- coefficients <- hollowCoefficients plan radius periodEnvelope coefficients'
      -- coefficients = getNormalizedHistogramArr $ hist
        -- normalizeFreqArr' std phiFreqs rhoFreqs .
  harmonicsArray <-
    createHarmonics
      numR2Freq
      phiFreq
      rhoFreq
      thetaFreq
      scaleFreq
      (-s)
      periodR2
      periodEnvelope
      coefficients
      -- harmonicsArray1 =
      --   createHarmonics
      --     numR2Freq
      --     phiFreq
      --     rhoFreq
      --     thetaFreq
      --     scaleFreq
      --     (-1.25)
      --     periodR2
      --     periodEnvelope
      --     coefficients
  let !initialSourceDistribution =
        computeInitialDistributionFourierPinwheel
          numR2Freq
          periodR2
          periodEnvelope
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          initSource
      !initialSinkDistribution =
        computeInitialDistributionFourierPinwheel
          numR2Freq
          periodR2
          periodEnvelope
          phiFreq
          rhoFreq
          thetaFreq
          scaleFreq
          initSink
      -- harmonicsArray' =
      --   centerHollow numR2Freq $
      --   pinwheelFourierCoefficientsAnatical
      --     numR2Freq
      --     phiFreq
      --     rhoFreq
      --     thetaFreq
      --     scaleFreq
      --     (-s)
      --     periodR2
      --     periodEnvelope
      -- harmonicsArray1' =
      --   -- L.map (centerHollowVector numR2Freq) $
      --   pinwheelFourierCoefficientsAnaticalList1
      --     numR2Freq
      --     phiFreq
      --     rhoFreq
      --     thetaFreq
      --     scaleFreq
      --     (-s)
      --     periodR2
      --     periodEnvelope
      -- harmonicsArray2' =
      --   -- L.map (centerHollowVector numR2Freq) $
      --   pinwheelFourierCoefficientsAnaticalList2
      --     numR2Freq
      --     phiFreq
      --     rhoFreq
      --     thetaFreq
      --     scaleFreq
      --     (-s)
      --     periodR2
      --     periodEnvelope
      -- initialSourceDistribution =
      --   computeInitialDistributionFull'
      --     numR2Freq
      --     periodR2
      --     phiFreq
      --     rhoFreq
      --     initSource
      -- initialSourceDistribution =
      --   computeInitialDistributionPowerMethodPinwheelBasis'
      --     numR2Freq
      --     (-s)
      --     periodR2
      --     phiFreq
      --     rhoFreq
      --     initSource
      -- initialSinkDistribution =
      --   computeInitialDistributionFull'
      --     numR2Freq
      --     periodR2
      --     phiFreq
      --     rhoFreq
      --     initSink
  -- source <-
  --   convoluvePinhweelBasisTest'
  --     coefficients
  --     harmonicsArray1'
  --     harmonicsArray2'
  --     initialSourceDistribution
  -- sink <-
  --   timeReversal' <$>
  --   convoluvePinhweelBasisTest'
  --     coefficients
  --     harmonicsArray1'
  --     harmonicsArray2'
  --     initialSinkDistribution
  -- let source =
  --       convoluvePinhweelBasis'
  --         coefficients
  --         harmonicsArray'
  --         initialSourceDistribution
  --     sink =
  --       timeReversal' $
  --       convoluvePinhweelBasis'
  --         coefficients
  --         harmonicsArray'
  --         initialSinkDistribution
  --     sourceMat =
  --       A.transpose . A.use . fromDFTArray . getDFTArrayVector $ source
  --     sinkMat = A.transpose . A.use . fromDFTArray . getDFTArrayVector $ sink-
  source <- convolve harmonicsArray initialSourceDistribution
  -- source1 <- convolve harmonicsArray1 initialSourceDistribution
  let numLogPolarFreq = ((2 * scaleFreq + 1) * (2 * thetaFreq + 1))
      -- sourceMat = A.transpose . toMatrixAcc $ source
      -- sourceMat1 = A.transpose . toMatrixAcc $ source1
      -- sourceMat2 =
      --   A.transpose . toMatrixAcc $
      --   parZipWithFPArray (VS.zipWith (-)) source source1
  -- sourceR2 <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freq
  --     numPointsRecon
  --     numLogPolarFreq
  --     periodR2
  --     deltaRecon
  --     numBatchR2
  --     sourceMat
  -- plotImageRepa (folderPath </> "SourcePower.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
  --   VU.map sqrt . toUnboxed . sumS . R.map (\x -> (magnitude x) ** 2) . rotate3D $
  --   sourceR2
  sourceR2 <- plotFPArray plan (folderPath </> "SourcePower.png") source
  -- source1R2 <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freq
  --     numPointsRecon
  --     numLogPolarFreq
  --     periodR2
  --     deltaRecon
  --     numBatchR2
  --     sourceMat1
  -- plotImageRepa (folderPath </> "SourcePower1.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
  --   VU.map sqrt . toUnboxed . sumS . R.map (\x -> (magnitude x) ** 2) . rotate3D $
  --   source1R2
  -- source2R2 <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freq
  --     numPointsRecon
  --     numLogPolarFreq
  --     periodR2
  --     deltaRecon
  --     numBatchR2
  --     sourceMat2
  -- plotImageRepa (folderPath </> "SourcePower2.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
  --   VU.map sqrt . toUnboxed . sumS . R.map (\x -> (magnitude x) ** 2) . rotate3D $
  --   source2R2
  sink <- convolve harmonicsArray initialSinkDistribution
  let sinkMat = A.transpose . toMatrixAcc $ sink
  -- completion <- completionField' plan source (timeReversal' sink)
  -- let completionMat = createCuMat . getDFTArrayVector $ completion
  -- print completionMat
  sinkR2 <- plotFPArray plan (folderPath </> "SinkPower.png") sink
  -- sinkR2 <-
  --   computeFourierSeriesR2StreamAcc
  --     ptxs
  --     numR2Freq
  --     numPointsRecon
  --     numLogPolarFreq
  --     periodR2
  --     deltaRecon
  --     numBatchR2
  --     sinkMat
  -- plotImageRepa (folderPath </> "SinkPower.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
  --   VU.map sqrt . toUnboxed . sumS . R.map (\x -> (magnitude x) ** 2) . rotate3D $
  --   sinkR2
  -- completionR2 <-
  --   computeFourierSeriesR2Stream
  --     deviceIDs
  --     ptxs
  --     numR2Freq
  --     numPointsRecon
  --     periodR2
  --     deltaRecon
  --     numBatchR2
  --     [transposeCuMat completionMat]
  printCurrentTime "Start ploting.."
  -- plotDFTArrayPower
  --   (folderPath </> "CoefficientsPower.png")
  --   numR2Freq
  --   numR2Freq
  --   source
  -- plotImageRepaComplex (folderPath </> "Coefficients.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numR2Freq :. numR2Freq) .
  --   VU.convert . L.foldl1' (VS.zipWith (+)) . getDFTArrayVector $
  --   source
  -- plotImageRepaComplex (folderPath </> "Source.png") .
  --   ImageRepa 8 .
  --   fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
  --   toUnboxed . sumS . rotate3D $
  --   sourceR2
  completionR2 <-
    completionFieldRepa
      plan
      (R.reshape
         (Z :. (L.length scaleFreqs) :. (L.length thetaFreqs) :. numPointsRecon :.
          numPointsRecon) $
       sourceR2)
      (timeReversalRepa thetaFreqs .
       R.reshape
         (Z :. (L.length scaleFreqs) :. (L.length thetaFreqs) :. numPointsRecon :.
          numPointsRecon) $
       sinkR2)
  plotImageRepa (folderPath </> "CompletionPower.png") .
    ImageRepa 8 .
    fromUnboxed (Z :. (1 :: Int) :. numPointsRecon :. numPointsRecon) .
    VU.map sqrt .
    toUnboxed . sumS . sumS . R.map (\x -> (magnitude x) ** 2) . rotate4D2 $
    completionR2
  printCurrentTime "Done"
  -- sparse
  let harmonicsArraySparse =
        computeHarmonicsArraySparse
          numPointsRecon
          deltaRecon
          numPointsRecon
          deltaRecon
          phiFreqs
          rhoFreqs
          thetaFreqs
          scaleFreqs
          -- (2 * pi)
          (log (periodR2 / 4))
          64
          s
      init =
        computeInitialDistribution'
          numPointsRecon
          numPointsRecon
          phiFreqs
          rhoFreqs
          0
          initSource
  harmonicsArrayDFT <-
    fmap (listArray (bounds harmonicsArraySparse)) .
    dftExecuteBatchP
      plan
      (DFTPlanID DFT1DG [numPointsRecon, numPointsRecon] [0, 1]) .
    L.map (VS.convert . toUnboxed . computeUnboxedS . makeFilter2D) . IA.elems $
    harmonicsArraySparse
  sourceSparse <- convolve' Source plan coefficients harmonicsArrayDFT init
  plotDFTArrayPower
    (folderPath </> "SourcePower_Sparse.png")
    numPointsRecon
    numPointsRecon
    sourceSparse

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-# LANGUAGE FlexibleContexts, NoMonomorphismRestriction #-}
module NeuralNetwork where
import Control.Arrow hiding (app)
import Control.Monad
import Control.Monad.Trans.State.Lazy
import Data.List
import Numeric.LinearAlgebra.HMatrix
import System.Random
import qualified Data.Vector.Storable as V

-- ActivationFunction wraps a function with its derivative, along with a name (for printing)
data ActivationFunction a = AF (a -> a) (a -> a) String
instance Show (ActivationFunction a) where show (AF _ _ n) = "ActivationFunction " ++ n

tanhAF = AF tanh (\x -> 1 - (tanh x)^2) "tanh"
logisticAF = AF f (\x -> f x * (1 - f x)) "logistic" where f x = (1 / (1 + exp (-x)))

data NeuralNetwork a = NN {
    getWeightMatrices :: [Matrix a],
    getActivationFunction :: ActivationFunction a
    } deriving Show

applyNN :: Numeric a => NeuralNetwork a -> Vector a -> Vector a
applyNN (NN mats (AF theta _ _)) x = foldl' ((cmap theta .) . flip (app . tr) . V.cons 1) x mats

-- (forwardPropagate nn x) returns the values at all the intermediate layers, before and after the activation function is applied
forwardPropagate :: Numeric a => NeuralNetwork a -> Vector a -> ([Vector a], [Vector a])
forwardPropagate (NN mats (AF theta _ _)) x = first tail . unzip $ scanl aux (undefined, V.cons 1 x) mats where
    aux (_,z) w = let v = app (tr w) z in (v, V.cons 1 (cmap theta v))

-- (backPropagate nn x y) returns all the errors (in the same shape as the weights) of applying nn to x, with target value y
backPropagate :: (Num (Vector a), Numeric a) => NeuralNetwork a -> Vector a -> Vector a -> [Matrix a]
backPropagate nn@(NN mats (AF theta theta' _)) x y = gradient where
    deltaStep (ds, d) (v, w) = let d' = (cmap theta' v) * (V.tail $ app w d) in (d':ds, d')
    (v:vs, zs) = first reverse $ forwardPropagate nn x
    ws = reverse mats
    dLast = 2 * (cmap theta v - y)
    ds = fst $ foldl' deltaStep ([dLast], dLast) (zip vs ws)
    gradient = zipWith outer zs ds

batchUpdate alpha nn@(NN oldWeights _) dataset = let
    (gradient, size) = foldr (\(x, y) (gr, n) -> (zipWith (+) gr (backPropagate nn x y), n+1)) (cycle [scalar 0], 0) dataset
    step = map (*scalar (alpha / size)) gradient
    newWeights = zipWith (-) oldWeights step
    in nn {getWeightMatrices = newWeights}

stochasticGradientDescent dataset alpha nn = foldl' (\nn' pt -> batchUpdate alpha nn' [pt]) nn dataset

initializeMatrix mkEntry (m, n) = runState (fmap (matrix m) $ replicateM (m*n) (state mkEntry))

randomlyWeightedNetwork seed dims af = NN (evalState (mapM mkMat (zip dims (tail dims))) (mkStdGen seed)) af where
    mkMat (m, n) = state $ initializeMatrix (randomR (-1, 1)) (n,m+1)


text:
module IllusoryContourShapePinwheelBasis where

import           Control.Monad                  as M
import           Control.Monad.Parallel         as MP
import qualified Data.Array.Accelerate          as A
import           Data.Array.Accelerate.LLVM.PTX
import           Data.Array.Repa                as R
import           Data.Binary
import           Data.Complex
import           Data.List                      as L
import           Data.Vector.Storable           as VS
import           Data.Vector.Unboxed            as VU
import           FokkerPlanck
import           Foreign.CUDA.Driver            as CUDA
import           FourierMethod.BlockMatrixAcc
import           FourierMethod.FourierSeries2D
import           Image.IO
import           Pinwheel.FourierSeries2D
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.Array
import           Utils.Time
import           Utils.List
import Data.Array.IArray as IA
import FourierPinwheel
import Utils.Distribution
import FourierPinwheel.GaussianEnvelopePinwheel
import Filter.Utils

main = do
  args@(deviceIDsStr:numPointsStr:deltaStr:thresholdStr:numPointsReconStr:deltaReconStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:tauStr:numTrailsStr:deltaTStr:poissonWeightStr:numR2FreqStr:periodR2Str:phiFreqsStr:rhoFreqsStr:thetaFreqsStr:scaleFreqsStr:initDistStr:initScaleStr:histFilePath:histFilePathCorner:stdR2Str:stdThetaStr:stdRStr:numBatchR2Str:numBatchR2FreqsStr:numBatchOriStr:batchSizeStr:sStr:writeFlagStr:numIterationStr:shape2DStr:radiusStr:numThreadStr:_) <-
    getArgs
  let deviceIDs = read deviceIDsStr :: [Int]
      numPoints = read numPointsStr :: Int
      delta = read deltaStr :: Double
      threshold = read thresholdStr :: Double
      numPointsRecon = read numPointsReconStr :: Int
      deltaRecon = read deltaReconStr :: Double
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      tau = read tauStr :: Double
      numTrails = read numTrailsStr :: Int
      deltaT = read deltaTStr :: Double
      poissonWeight = read poissonWeightStr :: Double
      numR2Freq = read numR2FreqStr :: Int
      periodR2 = read periodR2Str :: Double
      phiFreq = read phiFreqsStr :: Int
      phiFreqs = L.map fromIntegral [-phiFreq .. phiFreq]
      rhoFreq = read rhoFreqsStr :: Int
      rhoFreqs = L.map fromIntegral [-rhoFreq .. rhoFreq]
      thetaFreq = read thetaFreqsStr :: Int
      thetaFreqs = L.map fromIntegral [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Int
      scaleFreqs = L.map fromIntegral [-scaleFreq .. scaleFreq]
      initScale = read initScaleStr :: Double
      initDist = read initDistStr :: [(Double, Double, Double, Double)]
      initPoints = L.map (\(x, y, t, s) -> Point x y t s) initDist
      initSource = [L.head initPoints]
      initSink = [L.last initPoints]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/IllusoryContourShapePinwheelBasis"
      stdR2 = read stdR2Str :: Double
      stdTheta = read stdThetaStr :: Double
      stdR = read stdRStr :: Double
      numBatchR2 = read numBatchR2Str :: Int
      numBatchR2Freqs = read numBatchR2FreqsStr :: Int
      numBatchOri = read numBatchOriStr :: Int
      batchSize = read batchSizeStr :: Int
      s = read sStr :: Double
      writeFlag = read writeFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      shape2Ds = read shape2DStr :: [Points Shape2D]
      radius = read radiusStr :: Double
      periodEnv = periodR2 ^ 2 / 4 --  * sqrt 2 --  ^ 2 * 2
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  flagCorner <- doesFileExist histFilePathCorner
  hist <-
    case (getShape . L.head $ shape2Ds) of
      Circle _ _ ->
        if flag
          then do
            printCurrentTime "read coefficients from file"
            decodeFile histFilePath
          else do
            printCurrentTime "Start computing coefficients..."
            initialise []
            devs <- M.mapM device deviceIDs
            ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
            ptxs <- M.mapM createTargetFromContext ctxs
            sampleCartesian
              histFilePath
              folderPath
              ptxs
              numPoints
              periodEnv
              delta
              numOrientation
              initScale
              thetaSigma
              tau
              threshold
              s
              phiFreq
              rhoFreq
              thetaFreq
              scaleFreq
              stdR2
      KoffkaCross _ _ ->
        if flagCorner
          then do
            printCurrentTime "read coefficients from file"
            decodeFile histFilePathCorner
          else do
            printCurrentTime "Start computing coefficients..."
            initialise []
            devs <- M.mapM device deviceIDs
            ctxs <- M.mapM (\dev -> CUDA.create dev []) devs
            ptxs <- M.mapM createTargetFromContext ctxs
            sampleCartesianCorner
              histFilePathCorner
              folderPath
              ptxs
              numPoints
              periodEnv
              delta
              numOrientation
              initScale
              thetaSigma
              tau
              threshold
              s
              phiFreq
              rhoFreq
              thetaFreq
              scaleFreq
              stdR2
              poissonWeight
  printCurrentTime "Done"
  printCurrentTime "Start Convloution.."
  plan <-
    makePlan
      folderPath
      emptyPlan
      numPointsRecon
      numPointsRecon
      numR2Freq
      (2 * thetaFreq + 1)
      (2 * scaleFreq + 1)
      (2 * phiFreq + 1)
      (2 * rhoFreq + 1)
  printCurrentTime "Compute DFT Plan done."
  let coefficients = getNormalizedHistogramArr hist
      asteriskGaussianVec =
        case getShape (L.head shape2Ds) of
          Circle _ _ ->
            gaussianPinwheel
              numR2Freq
              periodR2
              stdR2
              10
              thetaFreq
              scaleFreq
              periodEnv
              stdTheta
              stdR
          KoffkaCross _ _ ->
            gaussianPinwheel1
              numR2Freq
              periodR2
              stdR2
              10
              thetaFreq
              scaleFreq
              periodEnv
              stdTheta
              stdR
  harmonicsArray <-
    createHarmonics
      numR2Freq
      phiFreq
      rhoFreq
      thetaFreq
      scaleFreq
      (-s)
      periodR2
      periodEnv
      coefficients
  M.mapM_
    (\shape2D -> do
       let points =
             L.map (\(x, y) -> Point (x) (y) 0 1) .
             getShape2DIndexList' . makeShape2D $
             shape2D
       printCurrentTime (show points)
       (bias, dftBias) --Full
                          -- (dftBias, bias) <- -- Discrete
          <-
         case getShape shape2D of
           Circle _ _ ->
             computeBiasFourierPinwheelFull
               plan
               numR2Freq
               thetaFreq
               scaleFreq
               (-s)
               periodR2
               periodEnv
               radius
               stdTheta
               stdR
               stdR2
               asteriskGaussianVec
               points
           KoffkaCross _ _ ->
             computeBiasFourierPinwheelKoffkaCorss
               plan
               numR2Freq
               thetaFreq
               scaleFreq
               (-s)
               periodR2
               periodEnv
               radius
               stdTheta
               stdR
               stdR2
               asteriskGaussianVec
               points
       let initDist =
             computeInitialDistributionPowerMethodFourierPinwheelFull
               numR2Freq
               phiFreq
               rhoFreq
               thetaFreq
               scaleFreq
               bias -- dftBias
       computeContourFourierPinwheel
         plan
         folderPath
         writeFlag
         harmonicsArray
         dftBias -- bias
         numIteration
         numBatchR2
         numPointsRecon
         deltaRecon
         periodR2
         periodEnv
         initDist
         (show . getShape $ shape2D)
         deviceIDs)
    shape2Ds


text:
{-# LANGUAGE OverloadedStrings, DeriveGeneric #-}

module Data.Wavelets.Scaling where
import Prelude hiding (length,replicate,map,maximum,minimum)
-- import Data.Wavelets
import Numeric.Statistics
import Statistics.Function
import Statistics.Sample
import GHC.Generics
import Data.Vector.Storable -- hiding (map)
import Data.Packed.Matrix
-- import Linear

{-| Scale factors are designed to make it possible to use a scaled wavelet to reconstruct an approximation
    of the original data series. ordinary least squares is used so that translation is available for approximation,
    in addition to scaling.

|-}

data SeriesFactors = SeriesFactors {
      seriesMin :: Double ,
      seriesMax :: Double ,
      seriesMean :: Double ,
      seriesCount  :: Int } deriving (Show,Generic)

computeSeriesFactors :: Vector Double -> SeriesFactors
computeSeriesFactors v = SeriesFactors mn mx avg tot
    where (mn, mx,avg, tot) = minMaxAndAverage v

newtype OldSeriesFactors = OSF SeriesFactors  deriving (Show)
newtype NewSeriesFactors = NSF SeriesFactors  deriving (Show)
newtype NewSeriesMatrix = NSM (Matrix Double) deriving (Show)
newtype OldSeriesMatrix = OSM (Matrix Double) deriving (Show)

-- | x scaling matrix is B in Y = XB + e matrix formulation of ordinary least squares

newtype ScalingMatrix = SM (Matrix Double) deriving (Show)

type SeriesMatrix = Matrix Double

-- ScaleMatrix is the scaling matrix for the least squares problem
computeSeriesMatrix :: SeriesFactors -> SeriesMatrix
computeSeriesMatrix (SeriesFactors mn mx amean _ ) = (fromColumns vlst)
    where onesList = replicate 3 1
          nsf' = fromList [mn,mx,amean]
          vlst = [onesList,nsf']

computeNewSeriesMatrix :: NewSeriesFactors -> NewSeriesMatrix
computeNewSeriesMatrix (NSF sf) = NSM $ computeSeriesMatrix sf

computeSeriesResult :: SeriesFactors -> SeriesMatrix
computeSeriesResult (SeriesFactors mn mx avg _ ) = (fromColumns vlst)
    where nsf' = fromList [mn,mx,avg]
          vlst = [nsf']

computeOldSeriesMatrix :: OldSeriesFactors -> OldSeriesMatrix
computeOldSeriesMatrix (OSF sf) = OSM $ computeSeriesResult sf

-- | ordinary least squares estimation for the multivariate model
--   Y = X B + e        rows are observations, columns are elements
--   mean e = 0, cov e = kronecker s I

computeScalingMatrix :: NewSeriesFactors -> OldSeriesFactors -> ScalingMatrix
computeScalingMatrix nsf osf = let (OSM _Y) = computeOldSeriesMatrix osf
                                   (NSM _X) = computeNewSeriesMatrix nsf
                                   (_B,_ ,_ ) = ols _X _Y
                               in SM _B

-- | compute the average of the absolute value of a foldable type by setting it up like
-- foldl (\(a,t) b -> absAvgItrFcn a b (t+1) ) (0,0) foldableThing
absAvgItrFcn ::(Num a, Fractional a)=> a -> a -> a -> (a,a)
absAvgItrFcn a b t = avgtuple
   where
     diff = (a - abs(b))
     dv  = abs(diff/t)
     tot  = a + dv
     avgtuple = (tot,t)

-- | compute the maximum absolute value of a foldable type by setting it up like
-- foldl (\mx b -> absMaxItrFcn mx b ) 0 foldableThing
absMaxItrFcn :: (Num a, Ord a) => a -> a -> a
absMaxItrFcn a b = max a (abs b)

-- | compute the max min and average of a Vector in one pass
-- (avg, max, min, tot)
minMaxAndAverage :: Vector Double -> (Double, Double ,Double , Int )
minMaxAndAverage v = (mn,mx,av,l)
    where
      av = mean v
      (mn, mx) = minMax v
      l = length v

applyScalingMatrix :: ScalingMatrix -> Vector Double -> Vector Double
applyScalingMatrix (SM sm) v = map scaleFcn v
    where tCoef = sm @@> (0,0)
          sCoef = sm @@> (1,0)
          scaleFcn x = sCoef * x + tCoef


text:
{-# LANGUAGE FlexibleInstances, BangPatterns #-}
-----------------------------------------------------------------------------
-- |
-- Module     : Foreign.VMath.VNum
-- Copyright  : Copyright (c) 2010, Patrick Perry <patperry@gmail.com>
-- License    : BSD3
-- Maintainer : Patrick Perry <patperry@gmail.com>
-- Stability  : experimental
--
-- Vector Num operations.
--

module Foreign.VMath.VNum
    where

import Foreign( Storable, Ptr, with, )
import Foreign.Storable.Complex()
import Data.Complex( Complex(..) )

import Foreign.BLAS.Double( dcopy )
import Foreign.VMath.Double
import Foreign.VMath.Zomplex

-- | Types with vectorized 'Num' operations.
class (Storable a, Num a) => VNum a where
    vShift :: Int -> a -> Ptr a -> Ptr a -> IO ()
    vAdd :: Int -> Ptr a -> Ptr a -> Ptr a -> IO ()
    vSub :: Int -> Ptr a -> Ptr a -> Ptr a -> IO ()
    vAxpby :: Int -> a -> Ptr a -> a -> Ptr a -> Ptr a -> IO ()

    vScale :: Int -> a -> Ptr a -> Ptr a -> IO ()
    vMul :: Int -> Ptr a -> Ptr a -> Ptr a -> IO ()

    vConj :: Int -> Ptr a -> Ptr a -> IO ()
    vNeg :: Int -> Ptr a -> Ptr a -> IO ()
    vAbs :: Int -> Ptr a -> Ptr a -> IO ()
    vSgn :: Int -> Ptr a -> Ptr a -> IO ()

instance VNum Double where
    vShift = vdShift
    {-# INLINE vShift #-}
    vAdd = vdAdd
    {-# INLINE vAdd #-}
    vSub = vdSub
    {-# INLINE vSub #-}
    vAxpby = vdAxpby
    {-# INLINE vAxpby #-}
    vScale = vdScale
    {-# INLINE vScale #-}
    vMul = vdMul
    {-# INLINE vMul #-}
    vConj n src dst =
        with (toEnum n) $ \pn ->
        with 1 $ \pone ->
            dcopy pn src pone dst pone
    {-# INLINE vConj #-}
    vNeg = vdNeg
    {-# INLINE vNeg #-}
    vAbs = vdAbs
    {-# INLINE vAbs #-}
    vSgn = vdSgn
    {-# INLINE vSgn #-}

instance VNum (Complex Double) where
    vShift n alpha pX pZ = with alpha $ \pAlpha ->
        vzShift n pAlpha pX pZ
    {-# INLINE vShift #-}
    vAdd = vzAdd
    {-# INLINE vAdd #-}
    vSub = vzSub
    {-# INLINE vSub #-}
    vAxpby n alpha pX beta pY pZ =
        with alpha $ \pAlpha ->
        with beta $ \pBeta ->
            vzAxpby n pAlpha pX pBeta pY pZ
    {-# INLINE vAxpby #-}
    vScale n alpha pX pZ = with alpha $ \pAlpha ->
        vzScale n pAlpha pX pZ
    {-# INLINE vScale #-}
    vMul = vzMul
    {-# INLINE vMul #-}
    vConj = vzConj
    {-# INLINE vConj #-}
    vNeg = vzNeg
    {-# INLINE vNeg #-}
    vAbs = vzAbs
    {-# INLINE vAbs #-}
    vSgn = vzSgn
    {-# INLINE vSgn #-}

max_stars_repo_name: patperry/hs-linear-algebra

max_issues_repo_name: cartazio/hs-cblas

max_forks_repo_name: cartazio/hs-cblas


text:
import Data.Complex

type CD = Complex Double

quadraticRoots :: (CD, CD, CD) -> (CD, CD)
quadraticRoots (a, b, c) =
    if   realPart b > 0
    then ((2*c) / (-b - d), (-b - d) / (2*a))
    else ((-b + d) / (2*a), (2*c) / (-b + d))
  where d = sqrt $ b^2 - 4*a*c

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData


text:
module MachineLearning.RandomTest
(
  tests
)

where

import Test.Framework (testGroup)
import Test.Framework.Providers.HUnit
import Test.HUnit
import Test.HUnit.Approx
import Test.HUnit.Plus

import MachineLearning.Random

import Data.List (nub)
import Control.Monad (foldM_)
import qualified Data.Vector as V
import qualified Data.Vector.Storable as SV
import qualified Numeric.LinearAlgebra as LA
import qualified System.Random as Rnd
import qualified Control.Monad.Random as RndM

sampleTest = do
  gen <- Rnd.newStdGen
  foldM_ sampleTestIter gen [1..25]

sampleTestIter gen i =
  let xs = V.fromList [1..100]
      n = 10 + i
      (ys, gen') = sample gen n xs
  in do
    assertEqual "uniqness" (V.length xs) (length . nub $ V.toList xs)
    assertEqual "length" n (V.length ys)
    assertBool "maximum" $ (V.maximum ys) <= (V.maximum xs)
    assertBool "minimum" $ (V.minimum ys) >= (V.minimum xs)
    assertEqual "uniqness of elements" n (length . nub $ V.toList ys)
    return gen'

randomRListTest = mapM_ (randomRListTestIter ((-1000, 1000)::(Int, Int))) [10..30]

randomRListTestIter range@(lo, hi) len = do
  rndList <- RndM.evalRandIO (getRandomRListM len range)
  assertEqual "length" len (length rndList)
  assertBool "minimum" $ lo <= (minimum rndList)
  assertBool "maximum" $ hi >= (maximum rndList)

randomRVectorTest = mapM_ (randomRVectorTestIter ((-2, 2))) [10..30]

randomRVectorTestIter range@(lo, hi) len = do
  rndVector <- RndM.evalRandIO (getRandomRVectorM len range)
  assertEqual "length" len (SV.length rndVector)
  assertBool "minimum" $ lo <= (SV.minimum rndVector)
  assertBool "maximum" $ hi >= (SV.maximum rndVector)

randomRMatrixTest = mapM_ (randomRMatrixTestIter ((-2, 2))) $ zip [10..15] [12..17]

randomRMatrixTestIter range@(lo, hi) (rows, cols) = do
  rndMatrix <- RndM.evalRandIO (getRandomRMatrixM rows cols range)
  assertEqual "rows" rows (LA.rows rndMatrix)
  assertEqual "columns" cols (LA.cols rndMatrix)
  assertBool "minimum" $ lo <= (LA.minElement rndMatrix)
  assertBool "maximum" $ hi >= (LA.maxElement rndMatrix)

tests = [ testCase "sample" sampleTest
        , testCase "getRandomRList" randomRListTest
        , testCase "getRandomRVector" randomRVectorTest
        , testCase "getRandomRMatrix" randomRMatrixTest
        ]


text:
{-# LANGUAGE TypeFamilies #-}
-- | Tools for implementing (and debugging the use of) gradient descent schemes.
module HordeAd.Core.OptimizerTools
  ( updateWithGradient
  , gradientIsNil, minimumGradient, maximumGradient
  , ArgsAdam(..), defaultArgsAdam
  , StateAdam(..), initialStateAdam
  , liftMatrix43, updateWithGradientAdam
  ) where

import Prelude

import           Control.Monad.ST.Strict (runST)
import qualified Data.Array.DynamicS as OT
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Generic.Mutable as VM
import           HordeAd.Internal.OrthotopeOrphanInstances (liftVT2)
import           Numeric.LinearAlgebra (Element, Matrix, Numeric, Vector)
import qualified Numeric.LinearAlgebra as HM
import           Numeric.LinearAlgebra.Data (flatten)
import           Numeric.LinearAlgebra.Devel
  (MatrixOrder (..), liftMatrix, liftMatrix2, matrixFromVector, orderOf)

import HordeAd.Internal.Delta (Domains, isTensorDummy)

{-
60% of heap allocation in matrix- and vector-based MNIST
with simple gradient descent (no mini-batches) is performed
by @updateWithGradient.updateVector@ below

  let updateVector i r = i - HM.scale gamma r

due to allocating once in @scale@ and again in @-@
(and there seems to be one more allocation judging by the numbers).
Something like the following code would be needed to eliminate
one allocation, but it would requre the hmatrix maintainer to expose
internal modules.

import           Internal.Vectorized
  ( FunCodeSV (Scale), FunCodeVV (Sub), applyRaw, c_vectorMapValR
  , c_vectorZipR, createVector )

minusTimesGamma :: Storable r => r -> Vector r -> Vector r -> Vector r
minusTimesGamma gamma u v = unsafePerformIO $ do
  r <- createVector (dim0 u)
  pval <- newArray [gamma]
  (v `applyRaw` (r `applyRaw` id))
    (c_vectorMapValR (fromei Scale) pval)
    #| "minusTimesGamma1"
  free pval
  (u `applyRaw` (v `applyRaw` (r `applyRaw` id)))
    (c_vectorZipR (fromei Sub))
    #| "minusTimesGamma2"
  return r

BTW, a version with HM.Devel.zipVectorWith makes
the test twice slower and allocate twice more

  let updateVector = zipVectorWith (\i r -> i - gamma * r)

and a version with Vector.Storable makes the test thrice slower
and allocate thrice more

  let updateVector = V.zipWith (\i r -> i - gamma * r)

which is probably a bug in stream fusion which, additionally in this case,
can't fuse with anything and so can't pay for its overhead.

-}

updateWithGradient :: (Numeric r, Floating (Vector r))
                   => r -> Domains r -> Domains r -> Domains r
updateWithGradient gamma (params0, params1, params2, paramsX)
                         (gradient0, gradient1, gradient2, gradientX) =
  let updateVector i r = i - HM.scale gamma r
      !params0New = updateVector params0 gradient0
      update1 i r = if V.null r  -- eval didn't update it, would crash
                    then i
                    else updateVector i r
      !params1New = V.zipWith update1 params1 gradient1
      update2 i r = if HM.rows r <= 0  -- eval didn't update it, would crash
                    then i
                    else liftMatrix2 updateVector i r
      !params2New = V.zipWith update2 params2 gradient2
      updateX i r = if isTensorDummy r  -- eval didn't update it, would crash
                    then i
                    else liftVT2 updateVector i r
                      -- TODO: this is slow; add @liftArray2@ and use HM,
                      -- unless we move away from HM; similarly other OT calls
      !paramsXNew = V.zipWith updateX paramsX gradientX
  in (params0New, params1New, params2New, paramsXNew)

gradientIsNil :: (Eq r, Numeric r)
              => Domains r -> Bool
gradientIsNil (gradient0, gradient1, gradient2, gradientX) =
  V.all (== 0) gradient0
  && V.all V.null gradient1
  && V.all (\r -> HM.rows r <= 0) gradient2
  && V.all isTensorDummy gradientX

minimumGradient :: (Ord r, Numeric r)
                => Domains r -> r
minimumGradient (gradient0, gradient1, gradient2, gradientX) =
  min (if V.null gradient0 then 0 else HM.minElement gradient0)
      (min (if V.null gradient1 then 0
            else V.minimum (V.map HM.minElement gradient1))
           (min (if V.null gradient2 then 0
                 else V.minimum (V.map HM.minElement gradient2))
                (if V.null gradientX then 0
                 else V.minimum (V.map OT.minimumA gradientX))))

maximumGradient :: (Ord r, Numeric r)
                => Domains r -> r
maximumGradient (gradient0, gradient1, gradient2, gradientX) =
  max (if V.null gradient0 then 0 else HM.maxElement gradient0)
      (max (if V.null gradient1 then 0
            else V.maximum (V.map HM.maxElement gradient1))
           (max (if V.null gradient2 then 0
                 else V.maximum (V.map HM.maxElement gradient2))
                (if V.null gradientX then 0
                 else V.maximum (V.map OT.maximumA gradientX))))

data ArgsAdam r = ArgsAdam
  { alpha   :: r
  , betaOne :: r
  , betaTwo :: r
  , epsilon :: r
  }

-- The defaults taken from
-- https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/Adam
defaultArgsAdam :: Fractional r => ArgsAdam r
defaultArgsAdam = ArgsAdam
  { alpha = 0.001
  , betaOne = 0.9
  , betaTwo = 0.999
  , epsilon = 1e-7
  }

data StateAdam r = StateAdam
  { tAdam :: Int  -- iteration count
  , mAdam :: Domains r
  , vAdam :: Domains r
  }

-- The arguments are just sample params0, for dimensions.
zeroParameters :: Numeric r
               => Domains r -> Domains r
zeroParameters (params0, params1, params2, paramsX) =
  let zeroVector v = runST $ do
        vThawed <- V.thaw v
        VM.set vThawed 0
        V.unsafeFreeze vThawed
  in ( zeroVector params0
     , V.map zeroVector params1
     , V.map (liftMatrix zeroVector) params2
     , V.map (\a -> OT.constant (OT.shapeL a) 0) paramsX )  -- fast allright

initialStateAdam :: Numeric r
                 => Domains r -> StateAdam r
initialStateAdam parameters0 =
  let zeroP = zeroParameters parameters0
  in StateAdam
       { tAdam = 0
       , mAdam = zeroP
       , vAdam = zeroP
       }

-- | Application of a vector function on the flattened matrices elements.
liftMatrix43 :: ( Numeric a, Numeric b, Numeric c, Numeric d
                , Element x, Element y, Element z )
             => (Vector a -> Vector b -> Vector c -> Vector d
                 -> (Vector x, Vector y, Vector z))
             -> Matrix a -> Matrix b -> Matrix c -> Matrix d
             -> (Matrix x, Matrix y, Matrix z)
liftMatrix43 f m1 m2 m3 m4 =
  let sz@(r, c) = HM.size m1
      rowOrder = orderOf m1
        -- checking @m4@ (gradient) makes RNN LL test much faster and BB slower
        -- so this needs much more benchmarking and understading to tweak
  in if sz == HM.size m2 && sz == HM.size m3 && sz == HM.size m4
     then
       let (vx, vy, vz) = case rowOrder of
             RowMajor -> f (flatten m1) (flatten m2) (flatten m3) (flatten m4)
             ColumnMajor -> f (flatten (HM.tr' m1)) (flatten (HM.tr' m2))
                              (flatten (HM.tr' m3)) (flatten (HM.tr' m4))
       in ( matrixFromVector rowOrder r c vx
          , matrixFromVector rowOrder r c vy
          , matrixFromVector rowOrder r c vz
          )
     else error $ "nonconformant matrices in liftMatrix43: "
                  ++ show (HM.size m1, HM.size m2, HM.size m3, HM.size m4)

-- TOOD: make sure this is not worse that OT.zipWith3A when transposing
-- between each application or that we never encounter such situations
--
-- | Application of a vector function on the flattened arrays elements.
liftArray43 :: ( Numeric a, Numeric b, Numeric c, Numeric d
               , Numeric x, Numeric y, Numeric z )
            => (Vector a -> Vector b -> Vector c -> Vector d
                -> (Vector x, Vector y, Vector z))
            -> OT.Array a -> OT.Array b -> OT.Array c -> OT.Array d
            -> (OT.Array x, OT.Array y, OT.Array z)
liftArray43 f m1 m2 m3 m4 =
  let sz = OT.shapeL m1
  in if sz == OT.shapeL m2 && sz == OT.shapeL m3 && sz == OT.shapeL m4
     then
       let (vx, vy, vz) = f (OT.toVector m1) (OT.toVector m2)
                            (OT.toVector m3) (OT.toVector m4)
       in ( OT.fromVector sz vx
          , OT.fromVector sz vy
          , OT.fromVector sz vz
          )
     else error
          $ "nonconformant arrays in liftArray43: "
            ++ show (OT.shapeL m1, OT.shapeL m2, OT.shapeL m3, OT.shapeL m4)

updateWithGradientAdam
  :: forall r. (Numeric r, Floating r, Floating (Vector r))
  => ArgsAdam r -> StateAdam r -> Domains r -> Domains r
  -> (Domains r, StateAdam r)
updateWithGradientAdam ArgsAdam{..}
                       StateAdam{ tAdam
                                , mAdam = (mAdam0, mAdam1, mAdam2, mAdamX)
                                , vAdam = (vAdam0, vAdam1, vAdam2, vAdamX)
                                }
                       (params0, params1, params2, paramsX)
                       (gradient0, gradient1, gradient2, gradientX) =
  let tAdamNew = tAdam + 1
      oneMinusBeta1 = 1 - betaOne
      oneMinusBeta2 = 1 - betaTwo
      updateVector :: Vector r -> Vector r
                   -> Vector r -> Vector r
                   -> (Vector r, Vector r, Vector r)
      updateVector mA vA p g =
        let mANew = HM.scale betaOne mA + HM.scale oneMinusBeta1 g
            vANew = HM.scale betaTwo vA + HM.scale oneMinusBeta2 (g * g)
            alphat = alpha * sqrt (1 - betaTwo ^ tAdamNew)
                             / (1 - betaOne ^ tAdamNew)
        in ( mANew
           , vANew
           , p - HM.scale alphat mANew
                 / (sqrt vANew + HM.scalar epsilon) )  -- the @scalar@ is safe
                      -- @addConstant@ would be better, but it's not exposed
      (!mAdam0New, !vAdam0New, !params0New) =
        updateVector mAdam0 vAdam0 params0 gradient0
      update1 mA vA p g = if V.null g  -- eval didn't update it, would crash
                          then (mA, vA, p)
                          else updateVector mA vA p g
      (!mAdam1New, !vAdam1New, !params1New) =
        V.unzip3 $ V.zipWith4 update1 mAdam1 vAdam1 params1 gradient1
      update2 mA vA p g = if HM.rows g <= 0  -- eval didn't update it; crash
                          then (mA, vA, p)
                          else liftMatrix43 updateVector mA vA p g
      (!mAdam2New, !vAdam2New, !params2New) =
        V.unzip3 $ V.zipWith4 update2 mAdam2 vAdam2 params2 gradient2
      updateX mA vA p g = if isTensorDummy g  -- eval didn't update it
                          then (mA, vA, p)
                          else liftArray43 updateVector mA vA p g
      (!mAdamXNew, !vAdamXNew, !paramsXNew) =
        V.unzip3 $ V.zipWith4 updateX mAdamX vAdamX paramsX gradientX
  in ( (params0New, params1New, params2New, paramsXNew)
     , StateAdam { tAdam = tAdamNew
                 , mAdam = (mAdam0New, mAdam1New, mAdam2New, mAdamXNew)
                 , vAdam = (vAdam0New, vAdam1New, vAdam2New, vAdamXNew)
                 }
     )

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
module Math.IRT.Model.TwoPLM
  ( TwoPLM (..)
  ) where

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.Model.FourPLM ( FourPLM(..) )
import Math.IRT.Model.Generic

data TwoPLM = TwoPLM { discrimination :: !Double
                     , difficulty     :: !Double
                     } deriving (Show)

instance Distribution TwoPLM where
    cumulative = cumulative . toFourPLM

instance ContDistr TwoPLM where
    density    = density . toFourPLM
    quantile _ = error "This shouldn't be needed"

instance DensityDeriv TwoPLM where
    densityDeriv = densityDeriv . toFourPLM

instance GenericModel TwoPLM where
    fromRasch         b = TwoPLM 1.0 b
    fromOnePLM        b = TwoPLM 1.7 b
    fromTwoPLM          = TwoPLM
    fromThreePLM  a b _ = TwoPLM   a b
    fromFourPLM a b _ _ = TwoPLM   a b

instance LogLikelihood TwoPLM where
    logLikelihood b = logLikelihood b . toFourPLM

toFourPLM :: TwoPLM -> FourPLM
toFourPLM (TwoPLM sa sb) = FourPLM sa sb 0.0 1.0

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE BangPatterns #-}
module Main where

import Criterion.Main
import Data.Massiv.Array as A
--import Data.Massiv.Array.SIMD
--import Data.Massiv.Array.Manifest.Vector as A
import Data.Massiv.Array.Unsafe as A
import Data.Massiv.Bench as A
import Prelude as P hiding ((<>))
import Data.Semigroup

--import Statistics.Matrix as S
--import Statistics.Matrix.Fast as SF

multArrsAlt :: Array U Ix2 Double -> Array U Ix2 Double -> Array U Ix2 Double
multArrsAlt arr1 arr2
  | n1 /= m2 =
    error $
    "(|*|): Inner array dimensions must agree, but received: " ++
    show (size arr1) ++ " and " ++ show (size arr2)
  | otherwise = compute $
    makeArrayR D (getComp arr1 <> getComp arr2) (Sz (m1 :. n2)) $ \(i :. j) ->
      A.foldlS (+) 0 (A.zipWith (*) (unsafeOuterSlice arr1 i) (unsafeOuterSlice arr2' j))
  where
    Sz2 m1 n1 = size arr1
    Sz2 m2 n2 = size arr2
    arr2' = computeAs U $ A.transpose arr2
{-# INLINE multArrsAlt #-}

main :: IO ()
main = do
  let !sz@(Sz2 _m _n) = Sz2 600 600
      !arr = arrRLightIx2 U Seq sz
      -- !arr2 = computeAs U arr
      -- !mat = S.Matrix m n $ A.toVector arr
      -- !mat' = S.transpose mat
      -- !arrV = arrRLightIx2 V Seq sz
      -- !arrV2 = computeAs V arrV
  defaultMain
    [ env (return (computeAs U (A.transpose arr))) $ \ arr' ->
        bgroup
          "Mult"
          [ bgroup
              "Seq"
              [ bench "(|*|)" $ whnfIO (setComp Seq arr |*| arr')
              -- , bench "multiplyTranspose" $
              --   whnf (computeAs U . multiplyTransposed (setComp Seq arr)) arr2
              -- , bench "multiplyTransposeSIMD" $
              --   whnf (computeAs P . multiplyTransposedSIMD arrV) arrV2
              , bench "multArrsAlt" $ whnf (multArrsAlt (setComp Seq arr)) arr'
              --, bench "multiply (dense-linear-algebra)" $ whnf (SF.multiply mat) mat'
              ]
          , bgroup
              "Par"
              [ bench "(|*|)" $ whnfIO (setComp Par arr |*| arr')
              , bench "fused (|*|)" $ whnfIO (setComp Par arr |*| A.transpose arr)
              -- , bench "multiplyTranspose" $
              --   whnf (computeAs U . multiplyTransposed (setComp Par arr)) arr2
              -- , bench "multiplyTransposeSIMD" $
              --   whnf (computeAs P . multiplyTransposedSIMD (setComp Par arrV)) arrV2
              , bench "multArrsAlt" $ whnf (multArrsAlt (setComp Par arr)) arr'
              ]
          ]
    ]


text:
{-|
Module      : GMM
Description : Definition of the Gaussian Mixture Model and Metropolis hatings sampler.
Copyright   : (c) Julian Kopka Larsen, 2015
Stability   : experimental

-}
module GMM (
    X,
    getElement
) where

import GHC.Arr (range)
import Data.List (foldl')
import Data.Array
import System.Random
import Numeric.LinearAlgebra
import Numeric.LinearAlgebra.Util hiding ((!))
import Debug.Trace

import Distributions -- (lnormalInvWishart, lMixNormWish, delta, Expr)
import Partition (Partition, Component, Move, naivefromNk, genMoves, applyMove, group)

import Math (X)

type Chain = [Partition]

type Likelihood = (X -> Double)

type Proposal = (Move, Double)

props :: Int -> Int -> Int -> [Proposal]
props seed n k = zip moves rand_accepts
                where rand_accepts = randomRs (0.0,1.0) (mkStdGen seed) :: [Double]
                      moves = genMoves seed n k

sample :: (Partition -> Move -> Double) -> Partition -> Proposal -> Partition
sample f prev_state (m, accept_prop)
    | accept_prop < f prev_state m = new_state -- trace ("(" ++ (show $ a!(node m)) ++ "->" ++ (show $ comp m) ++   ")" ++ "Accepted: " ++ show (f a m)) b
    | otherwise                    = prev_state -- trace ("(" ++ (show $ a!(node m)) ++ "->" ++ (show $ comp m) ++   ")" ++ "Rejected: " ++ show (f a m)) a
                                    where new_state = applyMove prev_state m

getElement :: X -> Int -> Int -> Int -> Partition
getElement x seed k nSamples = foldl' sampler start $ take nSamples (props seed n k)
                      where n = length x
                            start = naivefromNk n k
                            sampler = sample $ acceptanceRatio dq
                            dq = dlNormW x

acceptanceRatio :: (Partition -> Partition -> Double) -> Partition -> Move -> Double
acceptanceRatio dq x m = exp (dq x' x)
                        where x' = applyMove x m

max_stars_repo_name: juliankopkalarsen/FpStats

max_issues_repo_name: juliankopkalarsen/FpStats

max_forks_repo_name: juliankopkalarsen/FpStats


text:
{-# LANGUAGE QuasiQuotes, TemplateHaskell, ExplicitForAll, ScopedTypeVariables #-}
module AlgebraicDataTypes where

import Language.Rust.Inline
import Language.Rust.Inline.TH
import Language.Rust.Quote
import Test.Hspec

import Foreign.Marshal.Array
import Foreign.Ptr
import Foreign.Storable

import Data.Bifunctor ( Bifunctor(..) )
import Data.Char      ( toUpper )
import Data.Complex   ( Complex(..) )
import Data.Foldable  ( for_ )
import Data.Int       ( Int8, Int16, Int32, Int64 )

-- | A struct-like ADT where the fields have different sizes
data StructLike = StructLike  Int16 Int64 deriving (Show, Eq)
mkStorable [t| Storable StructLike |]

-- | A struct-like newtype ADT where the field is compound
newtype StructLike2 = StructLike2 (Int16, Int64) deriving (Show, Eq)
mkStorable [t| Storable StructLike2 |]

-- | An ADT where:
--
--   * Different variants have different alignment
--   * Variants have other complex structures
--
data Foo
  = Bar
  | Baz Char Int
  | Qux (Complex Float) Char
  deriving (Show, Eq)
mkStorable [t| Storable Foo |]

-- | An ADT where fields are nested ADTs
data Croc
  = Lob (Maybe Foo) Int
  | Boo Int8 Int8
  deriving (Show, Eq)
mkStorable [t| Storable Croc |]

-- | A polymorphic ADT. (From the @these@ package).
data These a b
  = This a
  | That b
  | Both a b
  deriving (Show, Eq)
mkStorable [t| forall a b. (Storable a, Storable b) => Storable (These a b) |]

-- | An ADT that needs more that a 'Word8' to store the tag
data Big a
  = C000 | C001 | C002 | C003 | C004 | C005 | C006 | C007 | C008 | C009
  | C010 | C011 | C012 | C013 | C014 | C015 | C016 | C017 | C018 | C019
  | C020 | C021 | C022 | C023 | C024 | C025 | C026 | C027 | C028 | C029
  | C030 | C031 | C032 | C033 | C034 | C035 | C036 | C037 | C038 | C039
  | C040 | C041 | C042 | C043 | C044 | C045 | C046 | C047 | C048 | C049
  | C050 | C051 | C052 | C053 | C054 | C055 | C056 | C057 | C058 | C059
  | C060 | C061 | C062 | C063 | C064 | C065 | C066 | C067 | C068 | C069
  | C070 | C071 | C072 | C073 | C074 | C075 | C076 | C077 | C078 | C079
  | C080 | C081 | C082 | C083 | C084 | C085 | C086 | C087 | C088 | C089
  | C090 | C091 | C092 | C093 | C094 | C095 | C096 | C097 | C098 | C099
  | C100 | C101 | C102 | C103 | C104 | C105 | C106 | C107 | C108 | C109
  | C110 | C111 | C112 | C113 | C114 | C115 | C116 | C117 | C118 | C119
  | C120 | C121 | C122 | C123 | C124 | C125 | C126 | C127 | C128 | C129
  | C130 | C131 | C132 | C133 | C134 | C135 | C136 | C137 | C138 | C139
  | C140 | C141 | C142 | C143 | C144 | C145 | C146 | C147 | C148 | C149
  | C150 | C151 | C152 | C153 | C154 | C155 | C156 | C157 | C158 | C159
  | C160 | C161 | C162 | C163 | C164 | C165 | C166 | C167 | C168 | C169
  | C170 | C171 | C172 | C173 | C174 | C175 | C176 | C177 | C178 | C179
  | C180 | C181 | C182 | C183 | C184 | C185 | C186 | C187 | C188 | C189
  | C190 | C191 | C192 | C193 | C194 | C195 | C196 | C197 | C198 | C199
  | C200 | C201 | C202 | C203 | C204 | C205 | C206 | C207 | C208 | C209
  | C210 | C211 | C212 | C213 | C214 | C215 | C216 | C217 | C218 | C219
  | C220 | C221 | C222 | C223 | C224 | C225 | C226 | C227 | C228 | C229
  | C230 | C231 | C232 | C233 | C234 | C235 | C236 | C237 | C238 | C239
  | C240 | C241 | C242 | C243 | C244 | C245 | C246 | C247 | C248 | C249
  | C250 | C251 | C252 | C253 | C254 | C255 | C256 | C257 | C258 | C259
  | C260 | C261 | C262 | C263 | C264 | C265 | C266 | C267 | C268 | C269
  | C270 | C271 | C272 | C273 | C274 | C275 | C276 | C277 | C278 | C279
  | C280 | C281 | C282 | C283 | C284 | C285 | C286 | C287 | C288 | C289
  | C290 | C291 | C292 | C293 | C294 | C295 | C296 | C297 | C298 | C299 a
  deriving (Show, Eq)
mkStorable [t| forall a. Storable a => Storable (Big a) |]

-- | An ADT with a mixture of polymorphism and not.
data Foo2 a
  = Bar2
  | Baz2 Word
  | Qux2 a a
  | Quux2 Int a
  deriving (Show, Eq)
mkStorable [t| forall a. Storable a => Storable (Foo2 a) |]

-- Set the context
extendContext basic
extendContext prelude
extendContext (singleton [ty| Cpx<f32> |] [t| Complex Float |])
extendContext (rustTyCtx [t| Foo |])
extendContext (rustTyCtx [t| StructLike |])
extendContext (rustTyCtx [t| StructLike2 |])
extendContext (rustTyCtx [t| Croc |])
extendContext (rustTyCtx [t| forall a b. These a b |])
extendContext (rustTyCtx [t| forall a. Big a |])
extendContext (rustTyCtx [t| forall a. Foo2 a |])
setCrateModule

[rust|
#[derive(Copy,Clone)]
#[repr(C)]
pub struct Cpx<T> { re: T, im: T }

use std::ops::Add;
impl<T: Add<Output=T>> Add for Cpx<T> {
  type Output = Cpx<T>;

  fn add(self, rhs: Self) -> Self {
    Cpx {
      re: self.re + rhs.re,
      im: self.im + rhs.im,
    }
  }
}
|]

-- | Function on 'Foo'
quux :: Foo -> Foo
quux (Baz 'a' 0) = Bar
quux (Baz c i) = Baz 'l' (i+1)
quux (Qux (_ :+ im) c) = Qux (9.0 :+ im) c
quux f = f

[rust|
impl Foo {
  pub fn quux(self) -> Self {
    match self {
      Foo::Baz('a', 0) => Foo::Bar,
      Foo::Baz(c, i) => Foo::Baz('l', i+1),
      Foo::Qux(Cpx { im, .. }, c) => Foo::Qux(Cpx { re: 9f32, im }, c),
      f => f,
    }
  }
}
|]

-- | Function on 'Croc'
croc :: Croc -> Croc
croc (Lob mfoo i) = Lob (fmap quux mfoo) (i + 4)
croc (Boo i1 i2) = Boo (negate i1) (abs i2)

[rust|
impl Croc {
  pub fn croc(self) -> Self {
    match self {
      Croc::Lob(mfoo, i) => Croc::Lob(mfoo.map(|x| x.quux()), i + 4),
      Croc::Boo(i1, i2) => Croc::Boo(-i1, i2.abs()),
    }
  }
}
|]

-- | Function 1 on 'StructLike'/'StructLike2'
in2 :: StructLike -> StructLike2
in2 (StructLike f i) = StructLike2 (f * 2, i * 3 - 8)

-- | Function 2 on 'StructLike'/'StructLike2'
out2 :: StructLike2 -> StructLike
out2 (StructLike2 (f, i)) = StructLike (f `quot` 2) ((i + 8) `quot` 3)

[rust|
impl StructLike {
  pub fn in2(self) -> StructLike2 {
    let StructLike(f,i) = self;
    StructLike2((f * 2, i * 3 - 8))
  }
}

impl StructLike2 {
  pub fn out2(self) -> StructLike {
    let StructLike2((f,i)) = self;
    StructLike(f / 2, (i + 8) / 3)
  }
}
|]

instance Bifunctor These where
  bimap f _ (This x) = This (f x)
  bimap _ g (That y) = That (g y)
  bimap f g (Both x y) = Both (f x) (g y)

[rust|
impl<A, B> These<A, B> {
  pub fn bimap<C, D, F, G>(self, f: F, g: G) -> These<C, D>
  where
    F: FnOnce(A) -> C,
    G: FnOnce(B) -> D,
  {
    match self {
      These::This(x) => These::This(f(x)),
      These::That(y) => These::That(g(y)),
      These::Both(x, y) => These::Both(f(x), g(y)),
    }
  }
}
|]

algebraicDataTypes :: Spec
algebraicDataTypes = describe "Algebraic data types" $ do
  it "Can marshal a `Complex Float` argument/return" $ do
    let z1, z2 :: Complex Float
        z1 = 1.3 :+ 4.5
        z2 = 6.7 :+ 8.9
    [rust| Cpx<f32> { $(z1: Cpx<f32>) + $(z2: Cpx<f32>) } |] `shouldBe` z1 + z2

  it "Can marshal a custom single-constructor ADT argument/return" $ do
    let s1 = StructLike 78 (negate 267)
        s2 = StructLike 92 45223
        s3 = StructLike2 (34, -92391)
        s4 = StructLike2 (576, 1234)

    for_ [s1,s2] $ \si ->
        [rust| StructLike2 { $(si: StructLike).in2() } |] `shouldBe` in2 si
    for_ [s3,s4] $ \si ->
        [rust| StructLike { $(si: StructLike2).out2() } |] `shouldBe` out2 si

  it "Can marshal a custom monomorphic ADT argument/return" $ do
    let f1, f2, f3, f4 :: Foo
        f1 = Baz 'a' 0
        f2 = Baz 'b' 2
        f3 = Qux (7.1 :+ 3.4) 'f'
        f4 = Bar

    for_ [f1,f2,f3,f4] $ \fi ->
      [rust| Foo { $(fi: Foo).quux() } |] `shouldBe` quux fi

  it "Can marshal nested monomorphic ADT arguments/returns" $ do
    let c1, c2, c3, c4, c5, c6, c7 :: Croc
        c1 = Lob (Just (Baz 'a' 0)) 2
        c2 = Lob (Just (Baz 'b' 2)) 6
        c3 = Lob (Just (Qux (7.1 :+ 3.4) 'f')) 8
        c4 = Lob (Just Bar) 9
        c5 = Lob Nothing 3
        c6 = Boo 3 (-2)
        c7 = Boo (-4) 2

    for_ [c1,c2,c3,c4,c5,c6,c7] $ \ci ->
      [rust| Croc { $(ci: Croc).croc() } |] `shouldBe` croc ci

  it "Can marshal polymorphic ADT arguments/returns" $ do
    let t1, t2, t3 :: These Int8 Int64
        t1 = This maxBound
        t2 = That 432442
        t3 = Both (maxBound - 3) 879

    for_ [t1,t2,t3] $ \ti ->
      let v1 = [rust| These<i16,i64> {
                 $(ti: These<i8,i64>).bimap(|x| x as i16 * 2, |y| y + 2)
               } |]
          v2 = bimap (\x -> fromIntegral x * 2) (+ 2) ti
      in v1 `shouldBe` v2

  it "Can marshal nested polymorphic ADT arguments/returns" $ do
    let t1, t2, t3 :: These (Maybe Int) (These Int8 Int8)
        t1 = This (Just 6)
        t2 = This Nothing
        t3 = That (This 8)
        t4 = That (That 9)
        t5 = That (Both 1 2)
        t6 = Both (Just 3) (That 8)
        t7 = Both Nothing (Both 3 5)
        t8 = Both (Just 213) (Both 78 98)

    for_ [t1,t2,t3,t4,t5,t6,t7,t8] $ \ti ->
      let v1 = [rust| These<Option<isize>,These<i8,i8>> {
                 $(ti: These<Option<isize>,These<i8,i8>>).bimap(
                    |oi| oi.map(|i| i + 5),
                    |t| t.bimap(|i| i + 2, |j| j * 3),
                 )
               } |]
          v2 = bimap (fmap (+5)) (bimap (+2) (*3)) ti
      in v1 `shouldBe` v2

  it "Can marshal a big ADT whose tag needs more than a `Word8`" $ do
    let b1, b2, b3, b4 :: Big Int64
        b1 = C000
        b2 = C160
        b3 = C298
        b4 = C299 89

    for_ [b1,b2,b3,b4] $ \bi ->
      let v1 = [rust| Big<i64> {
                  match $(bi: Big<i64>) {
                    Big::C160 => Big::C161,
                    Big::C299(i) => Big::C299(i+1),
                    b => b,
                  }
               } |]
          v2 = case bi of
                 C160 -> C161
                 C299 i -> C299 (i + 1)
                 b -> b
      in v1 `shouldBe` v2

  it "Can marshal a custom `Foo2 Int` and `Foo2 (Foo2 Int)` return" $ do
    let f1, f2, f3, f4 :: Foo2 Int
        f1 = Bar2
        f2 = Baz2 3
        f3 = Qux2 (-1) 2
        f4 = Quux2 (-8) 3

    let fooed f = case f of
                    Qux2 x y -> Qux2 (Qux2 x y) (Qux2 y x)
                    Quux2 i x -> Quux2 (i + 1) (Qux2 x x)
                    Bar2 -> Bar2
                    Baz2 w -> Baz2 w

    let fooed' f = [rust| Foo2<Foo2<isize>> {
      match $(f: Foo2<isize>) {
        Foo2::Qux2(x,y) => Foo2::Qux2(Foo2::Qux2(x,y), Foo2::Qux2(y,x)),
        Foo2::Quux2(i,x) => Foo2::Quux2(i+1, Foo2::Qux2(x,x)),
        Foo2::Bar2 => Foo2::Bar2,
        Foo2::Baz2(w) => Foo2::Baz2(w),
      }
    } |]

    fooed f1 `shouldBe` fooed' f1
    fooed f2 `shouldBe` fooed' f2
    fooed f3 `shouldBe` fooed' f3
    fooed f4 `shouldBe` fooed' f4

max_stars_repo_name: hexresearch/inline-rust

max_issues_repo_name: hexresearch/inline-rust

max_forks_repo_name: hexresearch/inline-rust


text:
{-# LANGUAGE BangPatterns #-}

-- | 2D FFT functions
module SatOpt.FFT where

import           Control.Parallel.Strategies   (parMap, rdeepseq)
import           Numeric.GSL.Fourier           (fft, ifft)
import           Numeric.LinearAlgebra.HMatrix

half [] = []
half (x:y:xs) = x : half xs

clone [] = []
clone (x:xs) = x : x : clone xs

-- | Take the FFT of a matrix of complex numbers
fft2dM :: Matrix ℂ -> Matrix ℂ
fft2dM m = fromColumns $ clone $ clone $ inter0
  where
    !inter0 = mapFFT $ inter1
    !inter1 = half $ half $ toColumns $ fromRows $ clone $ clone $ inter2
    !inter2 = mapFFT $ inter3
    !inter3 = half $ half $ toRows m
    mapFFT = parMap rdeepseq fft

-- | Take the inverse FFT of a matrix of complex numbers
ifft2dM :: Matrix ℂ -> Matrix ℂ
ifft2dM m = fromColumns $ mapIFFT $ toColumns inter
  where
    !inter = fromRows $ mapIFFT $ toRows m
    mapIFFT = parMap rdeepseq ifft

max_stars_repo_name: taktoa/SatOpt

max_issues_repo_name: taktoa/SatOpt

max_forks_repo_name: taktoa/SatOpt


text:
module Main where

import HopfieldMat as H
import Numeric.LinearAlgebra
import Codec.Picture
import Data.Either
import Control.Applicative

pixelsAsList :: Image PixelRGB8 -> [Int]
pixelsAsList img@(Image w h _) = go <$> [0..w-1] <*> [0..h-1]
  where go y x = pixelToBinary $ pixelAt img x y
        -- black is 1, everything else is zero
        pixelToBinary (PixelRGB8 r _ _) =
          fromIntegral $ 1 - signum r

listToImage :: Int -> Int -> [Int] -> Image PixelRGB8
listToImage w h d =  generateImage go w h
  where
    -- this is a horrible idea performance wise. yolo
    go x y = let v = d !! (y * w + x)
                 pv = fromIntegral $ 255 * (1 - v)
              in PixelRGB8 pv pv pv

loadPixels :: FilePath -> IO (Vector R)
loadPixels path = do
  img <- readImage path
  let pat = fromRight [] $ pixelsAsList . convertRGB8 <$> img
  let vpat = fromList $ fromIntegral <$> pat
  return vpat

main :: IO ()
main = do
  clean <- loadPixels "data/test.png"
  noisy <- loadPixels "data/test_noisy.png"
  let ws = H.train clean (H.initialWeights (size clean))
  putStrLn . show $ H.energy ws clean
  putStrLn . show $ H.energy ws noisy
  let newImg = H.feed ws noisy
  savePngImage "data/output_test.png" . ImageRGB8
    . listToImage 25 25 $ round <$> toList newImg

max_stars_repo_name: mfandl/hml

max_issues_repo_name: mfandl/hml

max_forks_repo_name: mfandl/hml


text:
import Data.Array.CArray
import Data.Complex
import Data.List (tails)
import Math.FFT (dft, idft)

convolution :: (Num a) => [a] -> [a] -> [a]
convolution x = map (sum . zipWith (*) (reverse x)) . spread
  where
    spread = init . tails . (replicate (length x - 1) 0 ++)

convolutionFFT :: [Complex Double] -> [Complex Double] -> [Complex Double]
convolutionFFT x y = elems $ idft $ liftArray2 (*) (fft x) (fft y)
  where
    fft a = dft $ listArray (1, length a) a

main :: IO ()
main = do
  let x = [1, 2, 1, 2, 1]
      y = [2, 1, 2, 1, 2]
  print $ convolution x y
  print $ convolutionFFT x y

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
{-# LANGUAGE GADTs, DataKinds, PolyKinds #-}
{-# LANGUAGE FlexibleContexts #-}

module TensorDAG where

import Data.Vinyl
import Data.Vinyl.Functor
import Data.Singletons
import Numeric.LinearAlgebra (Numeric)

import TensorHMatrix
import DAGIO
import VarArgs

makeGrad1 :: (a -> b -> a) -> HList '[a] -> Identity b -> HList '[a]
makeGrad1 g (Identity a :& RNil) (Identity b) = Identity (g a b) :& RNil

makeGrad2 :: (a -> b -> c -> (a, b)) -> HList '[a, b] -> Identity c -> HList '[a, b]
makeGrad2 g (Identity a :& Identity b :& RNil) (Identity c) = Identity a' :& Identity b' :& RNil
  where (a', b') = g a b c

makeDot :: Numeric a => Node (Tensor '[n] a) -> Node (Tensor '[n] a) -> IO (Node a)
makeDot = makeNode (uncurry2 dot, makeGrad2 gradDot)

makeMV :: (SingI n, IntegralN n, Usable a) => Node (Tensor '[n, m] a) -> Node (Tensor '[m] a) -> IO (Node (Tensor '[n] a))
makeMV = makeNode (uncurry2 mv, makeGrad2 gradMV)

makeMM :: (SingI n, IntegralN n, SingI k, Usable a) => Node (Tensor '[n, m] a) -> Node (Tensor '[m, k] a) -> IO (Node (Tensor '[n, k] a))
makeMM = makeNode (uncurry2 mm, makeGrad2 gradMM)

makeSelect i = makeNode (uncurry1 $ select i, makeGrad1 $ gradSelect i)


text:
{- |
Module      : UI.Tile
Description : Creation and execution of viewer tiles.
-}
module UI.Tile ( Tile()
               , renderTile
               , cancelTile
               , tileRect
               , ifModified
               , withSynchedTileBuffer
               ) where

import Task.Block
import Task.Concurrent
import Data.Planar

import Data.Color

import Control.Concurrent
import Control.Concurrent.Async

import Data.Complex
import Control.Monad (void)
import Data.Word
import Foreign.ForeignPtr
import Foreign.Ptr

newtype ImagePoint = ImagePoint (Double,Double)

instance Planar ImagePoint where
    toCoords (ImagePoint (x,y)) = (x, y)
    fromCoords = ImagePoint

-- | A tile in the image viewer.
data Tile = Tile
    { imageRect        :: Rectangle ImagePoint
      -- ^ The region in view space described by the tile.
    , tileBuffer       :: Synchronizable (ForeignPtr Word8)
      -- ^ The buffer into which the tile will draw.
    , tileWorker       :: Async ()
      -- ^ The worker thread which is drawing this tile.
    , shouldRedrawTile :: MVar ()
      -- ^ A value which signals that the tile needs to be redrawn.
    }

cancelTile :: Tile -> IO ()
cancelTile = void . forkIO . cancel . tileWorker -- cancel, but get on with life

withSynchedTileBuffer :: Tile -> (Ptr Word8 -> IO b) -> IO b
withSynchedTileBuffer tile action = synchedWith (tileBuffer tile) (`withForeignPtr` action)

-- | Unpack the width and height of a tile.
tileRect :: Tile -> (Int, Int)
tileRect tile = (floor w, floor h)
    where (w, h) = dimensions $ imageRect tile

-- | Perform an action, but only if the tile needs to be redrawn.
ifModified :: Tile -> IO () -> IO ()
ifModified tile f = do
    redraw <- tryTakeMVar $ shouldRedrawTile tile
    case redraw of
        Nothing -> return ()
        Just _  -> f

-- | Construct a tile from a dynamical system, and begin drawing to it.
renderTile :: (Word32 -> Word32 -> Complex Double -> Complex Double -> Ptr Word8 -> IO ()) -- ([(Double, Double)] -> IO [Color]) -- ^ The rendering action
           -> (Int, Int)   -- ^ The height and width of this tile.
           -> Rectangle (Double, Double)
              -- ^ The region of the dynamical plane corresponding
              --   to this tile.
           -> IO Tile      -- ^ An action which allocates the tile and
                           --   forks a task which draws into it.

renderTile renderingAction (width, height) mRect = do
    putStrLn ("renderTile at w=" ++ show width ++ " h=" ++ show height)

    buf <- mallocForeignPtrBytes (3 * (width + 16) * (height + 16))

    -- Initial fill of the image
    withForeignPtr buf $ \ptr ->
      sequence_ [ pokeColor ptr index grey | index <- [0 .. width * height - 1] ]

    let iRect = rectangle (ImagePoint (0,0))
                          (ImagePoint (fromIntegral width, fromIntegral height))

    redraw     <- newMVar ()  -- used to request a redraw
    managedBuf <- synchronized buf

    let (mRectWidth, mRectHeight) = dimensions mRect

    worker <- async $ progressively fillBlock
                    $ Block { coordToModel = convertRect iRect mRect . fromCoords
                            , compute = renderingAction
                            , logSampleRate = 1
                            , blockBuffer = managedBuf
                            , x0 = 0
                            , y0 = 0
                            , deltaX = mRectWidth / fromIntegral width
                            , deltaY = negate (mRectHeight / fromIntegral height)
                            , xStride = width
                            , xSize = width
                            , ySize = height
                            , shouldRedraw = redraw
                            }
    link worker

    return Tile { imageRect = iRect
                , tileBuffer = managedBuf
                , tileWorker = worker
                , shouldRedrawTile = redraw
                }

max_stars_repo_name: matt-noonan/FractalStream

max_issues_repo_name: matt-noonan/FractalStream

max_forks_repo_name: matt-noonan/FractalStream


text:
{-# LANGUAGE AllowAmbiguousTypes                      #-}
{-# LANGUAGE DataKinds                                #-}
{-# LANGUAGE FlexibleContexts                         #-}
{-# LANGUAGE GADTs                                    #-}
{-# LANGUAGE PartialTypeSignatures                    #-}
{-# LANGUAGE ScopedTypeVariables                      #-}
{-# LANGUAGE TupleSections                            #-}
{-# LANGUAGE TypeApplications                         #-}
{-# LANGUAGE TypeOperators                            #-}
{-# OPTIONS_GHC -fno-warn-partial-type-signatures     #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fplugin GHC.TypeLits.Normalise       #-}

import           Backprop.Learn
import           Control.DeepSeq
import           Control.Exception
import           Control.Monad
import           Control.Monad.IO.Class
import           Control.Monad.Trans.Class
import           Control.Monad.Trans.State
import           Data.Char
import           Data.Conduit
import           Data.Default
import           Data.Foldable
import           Data.Proxy
import           Data.Time
import           Data.Type.Equality
import           Data.Type.Tuple
import           GHC.TypeNats
import           Numeric.LinearAlgebra.Static.Backprop
import           Numeric.LinearAlgebra.Static.Vector
import           Numeric.Opto
import           System.Environment
import           Text.Printf
import qualified Conduit                               as C
import qualified Data.Conduit.Combinators              as C
import qualified Data.Set                              as S
import qualified Data.Text                             as T
import qualified Data.Vector.Sized                     as SV
import qualified Data.Vector.Storable.Sized            as SVS
import qualified System.Random.MWC                     as MWC
import qualified System.Random.MWC.Distributions       as MWC

-- | TODO: replace with 'LModel'
charRNN
    :: forall n h1 h2. (KnownNat n, KnownNat h1, KnownNat h2)
    => LModel _ _ (R n) (R n)
charRNN = fca softMax
       #: dropout @h2 0.25
       #: lstm
       #: dropout @h1 0.25
       #: lstm
       #: nilLM

oneHotChar
    :: KnownNat n
    => S.Set Char
    -> Char
    -> R n
oneHotChar cs = oneHotR . fromIntegral . (`S.findIndex` cs)

main :: IO ()
main = MWC.withSystemRandom @IO $ \g -> do
    sourceFile:_  <- getArgs
    charMap <- S.fromList <$> readFile sourceFile
    SomeNat (Proxy :: Proxy n) <- pure $ someNatVal (fromIntegral (length charMap))
    SomeNat (Proxy :: Proxy n') <- pure $ someNatVal (fromIntegral (length charMap - 1))
    Just Refl <- pure $ sameNat (Proxy @(n' + 1)) (Proxy @n)

    printf "%d characters found.\n" (natVal (Proxy @n))

    let model0 = charRNN @n @100 @50
        model  = trainState . unrollFinal @(SV.Vector 15) $ model0

    p0 <- initParamNormal model 0.2 g

    let report n b = do
          liftIO $ printf "(Batch %d)\n" (b :: Int)
          t0 <- liftIO getCurrentTime
          C.drop (n - 1)
          mp <- mapM (liftIO . evaluate . force) =<< await
          t1 <- liftIO getCurrentTime
          case mp of
            Nothing -> liftIO $ putStrLn "Done!"
            Just p@(p' :# s') -> do
              chnk <- lift . state $ (,[])
              liftIO $ do
                printf "Trained on %d points in %s.\n"
                  (length chnk)
                  (show (t1 `diffUTCTime` t0))
                let trainScore = testModelAll maxIxTest model (TJust p) chnk
                printf "Training error:   %.3f%%\n" ((1 - trainScore) * 100)

                forM_ (take 15 chnk) $ \(x,y) -> do
                  let primed = primeModel model0 (TJust p') x (TJust s')
                  testOut <- fmap reverse . flip execStateT [] $
                      iterateModelM ( fmap (oneHotR . fromIntegral)
                                    . (>>= \r -> r <$ modify (r:))    -- trace
                                    . (`MWC.categorical` g)
                                    . SVS.fromSized
                                    . rVec
                                    )
                            100 model0 (TJust p') y primed
                  printf "%s|%s\n"
                    (sanitize . (`S.elemAt` charMap) . fromIntegral . maxIndexR <$> (toList x ++ [y]))
                    (sanitize . (`S.elemAt` charMap) <$> testOut)
              report n (b + 1)

    C.runResourceT . flip evalStateT []
        . runConduit
        $ forever ( C.sourceFile sourceFile
                 .| C.decodeUtf8
                 .| C.concatMap T.unpack
                 .| C.map (oneHotChar charMap)
                 .| leadings
                  )
       .| skipSampling 0.02 g
       .| C.iterM (modify . (:))
       .| optoConduit
            def
            p0
            (adam def (modelGradStoch crossEntropy noReg model g))
       .| report 2500 0
       .| C.sinkNull

sanitize :: Char -> Char
sanitize c | isPrint c = c
           | otherwise = '#'

max_stars_repo_name: mstksg/backprop-learn

max_issues_repo_name: mstksg/backprop-learn

max_forks_repo_name: mstksg/backprop-learn