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

import Util

import Text.ParserCombinators.Parsec
import Text.ParserCombinators.Parsec.Number

import Data.String.Utils

import Data.String
import Data.CSV
import Data.List

import Data.Ord

import Data.Int

import Statistics.Sample as S hiding(stdDev)
import Data.Vector(fromList)

import Numeric

import Data.Maybe
import Data.Either
import qualified Data.Map.Strict as M

import Debug.Trace

import System.Environment

data Speedup = Speedup {
    num :: String,
    name :: String,
    time :: Double,
    nCores :: String,
    speedup :: Double,
    stdDev :: Double,
    factor :: Double,
    factorStdDev :: Double,
    overhead :: Double,
    stdDevForOverhead :: Double,
    runtimeX :: Double,
    runtimeY :: Double
} deriving (Show)

textBFLaTeX :: String -> String
textBFLaTeX str = "\"\\textbf{" ++ str ++ "}\""

textITLaTeX :: String -> String
textITLaTeX str = "\"\\textit{" ++ str ++ "}\""

rndSignificant :: Int -> Double -> Double
rndSignificant cnt val
      | val == 0.0 = 0.0
      | otherwise = (fromIntegral shifted) / magnitude
        where d::Double = fromIntegral (ceiling (logBase 10 (if val < 0 then -val else val)))
              power::Int64 = (fromIntegral cnt) - (floor d)

              magnitude::Double = fromIntegral (10 ^ power)
              shifted::Int64 = (round (val * magnitude))

formatOverheadForLaTeX :: Double -> String
formatOverheadForLaTeX overhead
   | overhead < 0 = textBFLaTeX (formatPercent overhead)
   | overhead > 0 = textITLaTeX (formatPercent overhead)

formatPercent :: Double -> String
formatPercent x = ((showFFloat Nothing (rndPercent x) "") ++ "\\%")

rndPercent :: Double -> Double
rndPercent x = (rndSignificant 2 (x * 100))

main :: IO ()
main = do
    args <- getArgs
    if (length args < 3)
    then
        do
            putStrLn $ "usage: <program> file1 file2 output [calculateWorst]"
    else
        do
            let (file1:file2:output:rest) = args

                lineToSpeedup :: [String] -> Maybe Speedup
                lineToSpeedup (num:name:numStrings) | (length numStrings) == 4 =
                                                        let [timeString, nCores, speedupString, stdDevString] = numStrings
                                                            maybeDoubles = convDoubles [timeString, speedupString, stdDevString]
                                                        in
                                                            if (isNothing maybeDoubles)
                                                            then
                                                                Nothing
                                                            else
                                                                let (Just [time, speedup, stdDev]) = maybeDoubles
                                                                in
                                                                    Just $ Speedup {
                                                                        num = num,
                                                                        name = name,
                                                                        time = time,
                                                                        nCores = nCores,
                                                                        speedup = speedup,
                                                                        stdDev = stdDev,
                                                                        factor = 1,
                                                                        factorStdDev = 0,
                                                                        overhead = 0,
                                                                        stdDevForOverhead = 0,
                                                                        runtimeX = 0,
                                                                        runtimeY = 0
                                                                    }
                                                   | otherwise = Nothing
                lineToSpeedup _ = Nothing

                diff :: Speedup -> Speedup -> Maybe Speedup
                diff x y
                    | nCores x == nCores y && num x == num y = Just $ Speedup {
                            num = num x,
                            name = (name x) ++ "-" ++ (name y),
                            time = (time x) - (time y),
                            nCores = nCores x,
                            speedup = (speedup x) - (speedup y),
                            stdDev = max (stdDev x) (stdDev y),
                            -- care: different order than time. cause reasons
                            factor = (time y) / (time x),
                            factorStdDev = (max (stdDev x) (stdDev y)) / (time x),
                            overhead = (time y - time x) / (time y),
                            stdDevForOverhead = (max (stdDev x) (stdDev y)) / (time y),
                            runtimeX = time x,
                            runtimeY = time y
                        }
                    | otherwise = Nothing

                legend :: String
                legend = "\"\",\"name\",\"time\",\"nCores\",\"speedup\",\"max stddev\",\"factor\",\"factorStdDev\",\"overhead\",\"stdDevForOverhead\",\"runtimeX\",\"runtimeY\"" ++ "\n"

                str :: String -> String
                str st = "\"" ++ st ++ "\""

                toStringLn :: Speedup -> String
                toStringLn x = (toString x) ++ "\n"

                toString :: Speedup -> String
                toString x = (str (num x)) ++ "," ++ (str $ name x) ++ "," ++  (show $ time x) ++ "," ++
                                    (nCores x) ++ "," ++ (show $ speedup x) ++ "," ++ (show $ stdDev x) ++ "," ++
                                    (show $ factor x) ++ "," ++ (show $ factorStdDev x) ++ "," ++ (show $ overhead x) ++
                                    "," ++ (show $ stdDevForOverhead x) ++ "," ++ (show $ runtimeX x) ++  "," ++ (show $ runtimeY x)

                handleParse :: Either ParseError [[String]] -> [Speedup]
                handleParse (Right lines) = catMaybes $ map lineToSpeedup lines
                handleParse _  = []

            linesOrError1 <- parseFromFile csvFile file1
            linesOrError2 <- parseFromFile csvFile file2

            let speedups1 = handleParse $ linesOrError1
                speedups2 = handleParse $ linesOrError2
            if (length speedups1 == 0 || length speedups2 == 0 || length speedups1 /= length speedups2)
            then
                putStrLn "parse Error!"
            else
                do
                    let maybeDiffs = zipWith diff speedups1 speedups2
                        diffs = {-traceShowId $-} catMaybes maybeDiffs
                    if (length diffs /= length maybeDiffs)
                    then
                        putStrLn $ "parse Error!"
                    else
                        if (length args >= 4)
                        then
                            let
                                -- hacky, but nvm
                                (_:restOfRest) = rest

                                readBool :: String -> Bool
                                readBool = read

                                calcOutput :: Bool -> Speedup
                                calcOutput True = maximumBy (comparing factor) diffs
                                calcOutput False = minimumBy (comparing factor) diffs

                                rnd :: Speedup -> Speedup
                                rnd x = Speedup {
                                                    num = num x,
                                                    name = name x,
                                                    time = time x,
                                                    nCores = nCores x,
                                                    speedup = speedup x,
                                                    stdDev = stdDev x,
                                                    factor = rndVal (factor x),
                                                    factorStdDev = factorStdDev x,
                                                    overhead = rndVal $ overhead x,
                                                    stdDevForOverhead = stdDevForOverhead x,
                                                    runtimeX = runtimeX x,
                                                    runtimeY = runtimeY x
                                                }

                                rndVal = rndSignificant 3

                                meanMeanOverhead = S.mean . fromList

                                roundedGeometricMean = rndVal . S.geometricMean . fromList

                                roundedMean = rndVal . S.mean . fromList
                            in
                                do appendFile output $ (toString $ rnd $ calcOutput True)
                                    ++ ","
                                    ++ (toString $ rnd $ calcOutput False)
                                    ++ ","
                                    -- geometric mean of factor
                                    ++ (showFFloat Nothing $ roundedMean $ map (factor) diffs) ""
                                    ++ ","
                                    --
                                    ++ (showFFloat Nothing $ rndVal $ factorStdDev $ maximumBy (comparing factorStdDev) diffs) ""
                                    ++ ","
                                    ++ (formatOverheadForLaTeX $ meanMeanOverhead $ map (overhead) diffs)
                                    ++ ","
                                    ++ (formatPercent $ stdDevForOverhead $ maximumBy (comparing stdDevForOverhead) diffs)
                                    ++ ","
                                    ++ (showFFloat Nothing $ rndVal $ runtimeX $ maximumBy (comparing nCores) diffs) ""
                                    ++ "\n"

                        else
                            do writeFile output $ legend ++ (concat $ map toStringLn diffs)

max_stars_repo_name: Parrows/Parrows

max_issues_repo_name: Parrows/Parrows

max_forks_repo_name: Parrows/Parrows


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

module Geometry.Matrix
(
    -- * Functions
    getIndependent
)

where

import Numeric.LinearAlgebra.HMatrix
import qualified Data.Matrix as M
import Data.Function
-- instance Element Rational

-- getIndependent :: [[Rational]] -> [[Rational]] -> Int -> Int -> M.Matrix Rational-> [[Rational]]
-- getIndependent initial accum currRank maxRank slackMat
--     | currRank == maxRank = accum
--     | accum == [] = getIndependent (tail initial) (return $ head initial) 1 maxRank slackMat
--     | otherwise =   if newRank > currRank
--                     then getIndependent (tail initial) (candidate:accum) (currRank+1) maxRank slackMat
--                     else getIndependent (tail initial) accum currRank maxRank slackMat
--         where
--             slack = map (map fromRational) $ M.toLists slackMat
--             candidate = head initial
--             setVectors = map (1:) $ candidate:accum
--             matrix =  fromLists $ zipRev ((map (map fromRational)) setVectors :: [[Double]]) (slack)
--             newRank = rank matrix

-- | Get independet rows/columns of a non square matrix
getIndependent :: [[Rational]] -> [[Rational]] -> Int -> Int -> M.Matrix Rational-> [[Rational]]
getIndependent initial accum currRank maxRank slackMat
    | currRank == maxRank = accum
    | accum == [] = getIndependent (tail initial) (return $ head initial) 1 maxRank slackMat
    | otherwise =   if newRank > currRank
                    then getIndependent (tail initial) (candidate:accum) (currRank+1) maxRank slackMat
                    else getIndependent (tail initial ++ [candidate]) accum currRank maxRank slackMat
        where
            slack = map (map fromRational) $ M.toLists slackMat
            candidate = head initial
            setVectors = map (1:) $ candidate:accum
            matrix =  fromLists $ zipRev ((map (map fromRational)) setVectors :: [[Double]]) (slack)
            newRank = rank matrix

zipRev :: [[Double]] -> [[Double]] -> [[Double]]
zipRev a = (reverse . ((zipWith (++)) `on` reverse) a)

ext: hs

max_stars_repo_name: matheusvms/facet-enumeration

max_issues_repo_name: matheusvms/facet-enumeration

max_forks_repo_name: matheusvms/facet-enumeration


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

type Vector = CArray Int (Complex Double)

(.*), (.+) :: Vector -> Vector -> Vector
a .* b = liftArray2 (*) a b
a .+ b = liftArray2 (+) a b

normalize :: Double -> Vector -> Vector
normalize dx v =
  let factor = 1 / sqrt dx / norm2 v :+ 0
   in liftArray (factor *) v

data Parameters = Parameters
  { xmax :: Double
  , res :: Int
  , dt :: Double
  , timesteps :: Int
  , dx :: Double
  , x :: Vector
  , dk :: Double
  , ks :: Vector
  , imTime :: Bool
  }

defaultParameters :: Parameters
defaultParameters = makeParameters 10 512 0.01 1000 True

makeParameters :: Double -> Int -> Double -> Int -> Bool -> Parameters
makeParameters xmax res dt timesteps imTime =
  let fi = fromIntegral
      rng = (0, res - 1)
      ks = [0 .. div res 2 - 1] ++ [-div res 2 .. -1]
   in Parameters
        xmax
        res
        dt
        timesteps
        (2 * xmax / fi res)
        (listArray rng $
         map (\n -> xmax * (-1 + 2 * fi n / fi res) :+ 0) [1 .. res])
        (pi / xmax)
        (listArray rng $ map ((:+ 0) . (pi / xmax *) . fi) ks)
        imTime

data Operators = Operators
  { v :: Vector
  , rStep :: Vector
  , kStep :: Vector
  , wfc :: Vector
  }

makeOperators :: Parameters -> Complex Double -> Complex Double -> Operators
makeOperators param v0 wfc0 =
  let rng = (0, res param - 1)
      time
        | imTime param = dt param :+ 0
        | otherwise = 0 :+ dt param
      v = liftArray (\x -> 0.5 * (x - v0) ^ 2) (x param)
      rStep = liftArray (\x -> exp (-0.5 * time * x)) v
      kStep = liftArray (\k -> exp (-0.5 * time * k ^ 2)) (ks param)
      wfc = liftArray (\x -> exp (-(x - wfc0) ^ 2 / 2)) (x param)
   in Operators v rStep kStep (normalize (dx param) wfc)

evolve :: Parameters -> Operators -> [Operators]
evolve param op@(Operators _ rStep kStep _) = iterate splitop op
  where
    splitop op = op {wfc = wfc' op}
    wfc' = norm . (rStep .*) . idft . (kStep .*) . dft . (rStep .*) . wfc
    norm = if imTime param then normalize (dx param) else id

calculateEnergy :: Parameters -> Operators -> Double
calculateEnergy param ops = (* dx param) . sum . map realPart $ elems totalE
  where
    totalE = potentialE .+ kineticE
    potentialE = wfcConj .* v ops .* wfc ops
    kineticOp = liftArray ((/ 2) . (^ 2)) (ks param)
    kineticE = wfcConj .* idft (kineticOp .* dft (wfc ops))
    wfcConj = liftArray conjugate $ wfc ops

-- Use gnuplot to make an animated  GIF using ../gnuplot/plot_output.plt
-- $ gnuplot -e "folder='../haskell'" plot_output.plt
printEvolution :: Parameters -> [Operators] -> IO ()
printEvolution param =
  mapM_ (export . (format <$>)) . zip [0 ..] . take 100 . skip
  where
    skip (x:xs) = x : skip (drop (div (timesteps param) 100 - 1) xs)
    format (Operators v _ _ wfc) =
      let density = liftArray ((^ 2) . abs) wfc
          values = map (map (show . realPart) . elems) [x param, density, v]
       in intercalate "\n" $ map (intercalate "\t") $ transpose values
    export (i, f) = writeFile ("output" ++ pad (show i) ++ ".dat") f
    pad n = replicate (5 - length n) '0' ++ n

main :: IO ()
main = do
  let p = defaultParameters
      o = makeOperators p 0 4
      evol = evolve p o
  print $ calculateEnergy p (evol !! timesteps p)
  printEvolution p evol

max_stars_repo_name: atocil/algorithm-archive

max_issues_repo_name: atocil/algorithm-archive

max_forks_repo_name: atocil/algorithm-archive


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}

{-# LANGUAGE CPP                   #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE FlexibleContexts      #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE GADTs                 #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NoStarIsType          #-}
{-# LANGUAGE OverloadedLabels      #-}
{-# LANGUAGE OverloadedStrings     #-}
{-# LANGUAGE PolyKinds             #-}
{-# LANGUAGE RankNTypes            #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-# LANGUAGE UndecidableInstances  #-}

import           Control.Monad
import           System.Random.MWC

import           Data.Constraint              (Dict (..))
import           Data.List                    (zipWith5)
import           Data.Proxy
import           Data.Reflection
import           Data.Singletons
import           Data.Singletons.TypeLits     hiding (natVal)
import           GHC.TypeLits
import           Unsafe.Coerce                (unsafeCoerce)

import           Numeric.LinearAlgebra.Static (L, R)
import qualified Numeric.LinearAlgebra.Static as H

import           Criterion.Main

import           Grenade
import           Grenade.Utils.LinearAlgebra
import           Grenade.Utils.ListStore

main :: IO ()
main = do
  defaultMain
    [ bgroup
        "batchnorm forward (CxHxW)"
        [ benchBatchNorm "1x1x64" 1 1 64
        , benchBatchNorm "1x64x64" 1 64 64
        , benchBatchNorm "64x64x64" 64 64 64
        ]
    , bgroup
        "convolutions with bias (CxHxW)"
        [ benchBiasConvolution "1024x64x64, 125 kernels" 1024 64 64 125 1 1 64 64
        , benchBiasConvolution "3x416x416, 64 kernels" 3 416 416 64 2 1 208 208
        ]
    , bgroup
        "leaky relu (CxHxW)"
        [ benchLeakyRelu "3x416x416" 3 416 416
        , benchLeakyRelu "1024x16x16" 1024 16 16
        ]
    ]

-- BENCHMARK GENERATION FUNCTIONS

benchBiasConvolution :: String -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> Benchmark
benchBiasConvolution name channels width height filters strides kernels outWidth outHeight
  | strides * (outHeight - 1) <= (height - kernels + 1) - 1 && (height - kernels + 1) <= (outHeight * strides) && strides * (outWidth - 1) <= (width - kernels + 1) - 1 && (width - kernels + 1) <= (outWidth * strides)
      = case (someNatVal (fromIntegral channels), someNatVal (fromIntegral width), someNatVal (fromIntegral height), someNatVal (fromIntegral filters), someNatVal (fromIntegral strides), someNatVal (fromIntegral kernels), someNatVal (fromIntegral outWidth), someNatVal (fromIntegral outHeight)) of
          (Just (SomeNat (Proxy :: Proxy channels)), Just (SomeNat (Proxy :: Proxy width)), Just (SomeNat (Proxy :: Proxy height)), Just (SomeNat (Proxy :: Proxy filters)), Just (SomeNat (Proxy :: Proxy strides)), Just (SomeNat (Proxy :: Proxy kernels)), Just (SomeNat (Proxy :: Proxy outWidth)), Just (SomeNat (Proxy :: Proxy outHeight))) ->
            case (channels, width, height
                 , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strides * (outHeight - 1) <= (height - kernels) ) )
                 , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (height - kernels) <= (outHeight * strides) - 1 ) )
                 , (unsafeCoerce (Dict :: Dict ()) :: Dict ( strides * (outWidth - 1) <= (width - kernels ) ) )
                 , (unsafeCoerce (Dict :: Dict ()) :: Dict ( (width - kernels) <= (outWidth * strides) - 1 ) ) ) of
              (1, 1, _, _, _, _, _)
                -> error "1D convolutions are not allowed"
              (1, _, _, Dict, Dict, Dict, Dict)
                -> env (generateBiasConvEnv :: IO (Convolution 'WithBias 'NoPadding 1 filters kernels kernels strides strides, S ('D2 height width))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D2 height width) -> S ('D3 outHeight outWidth filters )) x
              (_, _, _, Dict, Dict, Dict, Dict)
                -> env (generateBiasConvEnv :: IO (Convolution 'WithBias 'NoPadding channels filters kernels kernels strides strides, S ('D3 height width channels))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D3 height width channels) -> S ('D3 outHeight outWidth filters)) x
  where
    generateBiasConvEnv :: forall channels filters kernel1 kernel2 strides1 strides2 s.
                            ( SingI s, KnownNat channels, KnownNat filters, KnownNat kernel1, KnownNat kernel2,
                              KnownNat strides1, KnownNat strides2 )
                            => IO (Convolution 'WithBias 'NoPadding channels filters kernel1 kernel2 strides1 strides2, S s)
    generateBiasConvEnv = do
      x     <- randomOfShape
      gen   <- createSystemRandom
      layer <- createRandomWith UniformInit gen
      return (layer, x)

benchBatchNorm :: String -> Int -> Int -> Int -> Benchmark
benchBatchNorm name channels width height
  = case (someNatVal (fromIntegral channels), someNatVal (fromIntegral width), someNatVal (fromIntegral height)) of
          (Just (SomeNat (Proxy :: Proxy channels)), Just (SomeNat (Proxy :: Proxy width)), Just (SomeNat (Proxy :: Proxy height))) ->
            case (channels, width, height) of
              (1, 1, _) -> env (generateBatchNormEnv False :: IO (BatchNorm 1 1 width 90, S ('D1 width))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D1 width) -> S ('D1 width)) x
              (1, _, _) -> env (generateBatchNormEnv False :: IO (BatchNorm 1 width height 90, S ('D2 width height))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D2 width height) -> S ('D2 width height)) x
              (_, _, _) -> env (generateBatchNormEnv False :: IO (BatchNorm channels width height 90, S ('D3 width height channels))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D3 width height channels) -> S ('D3 width height channels)) x
  where
    generateBatchNormEnv :: forall c h w s. (KnownNat c, KnownNat h, KnownNat w, SingI s)
                         => Bool -> IO (BatchNorm c h w 90, S s)
    generateBatchNormEnv training = do
      x     <- randomOfShape
      gens  <- replicateM 4 createSystemRandom
      seeds <- mapM uniform gens :: IO [Int]
      let [gamma, beta, running_mean, running_var] = map (\s -> H.randomVector s H.Uniform) seeds :: [R c]
          ε            = 0.00001
      return (BatchNorm training (BatchNormParams gamma beta) running_mean running_var ε mkListStore, x)

benchLeakyRelu :: String -> Int -> Int -> Int -> Benchmark
benchLeakyRelu name channels width height
  = case (someNatVal (fromIntegral channels), someNatVal (fromIntegral width), someNatVal (fromIntegral height)) of
          (Just (SomeNat (Proxy :: Proxy channels)), Just (SomeNat (Proxy :: Proxy width)), Just (SomeNat (Proxy :: Proxy height))) ->
            case (channels, width, height) of
              (1, 1, _) -> env (generateLeakyReluEnv :: IO (LeakyRelu, S ('D1 width))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D1 width) -> S ('D1 width)) x
              (1, _, _) -> env (generateLeakyReluEnv :: IO (LeakyRelu, S ('D2 width height))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D2 width height) -> S ('D2 width height)) x
              (_, _, _) -> env (generateLeakyReluEnv :: IO (LeakyRelu, S ('D3 width height channels))) $ \ ~(layer, x) -> bench name $ nf (snd . runForwards layer :: S ('D3 width height channels) -> S ('D3 width height channels)) x
  where
    generateLeakyReluEnv :: forall s. (SingI s) => IO (LeakyRelu, S s)
    generateLeakyReluEnv = do
      x     <- randomOfShape
      gen   <- createSystemRandom
      alpha <- uniformR (-1, 1) gen :: IO RealNum
      return (LeakyRelu alpha, x)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE TemplateHaskell #-}
{-# LANGUAGE DataKinds #-}
module Lib
    (
    ) where
import Control.Lens
import Numeric.LinearAlgebra.Static as NLS
import Numeric.Backprop
import Debug.SimpleReflect as SR
data Network = Net { _weight1 :: NLS.L 20 100
                   , _bias1   :: NLS.R 20
                   , _weight2 :: NLS.L  5  20
                   , _bias2   :: NLS.R  5
                   }
Control.Lens.makeLenses ''Network
runNet :: Network -> NLS.R 100 -> NLS.R 5
runNet net x = z
  where
    y = logistic $ (net ^. weight1) #> x + (net ^. bias1)
    z = logistic $ (net ^. weight2) #> y + (net ^. bias2)
logistic :: Floating a => a -> a
logistic x = 1 / (1 + exp (-x))
myFun :: Floating a => a -> a -> a
myFun x y = x + y
myFunc :: Floating a => a -> a
myFunc x = sqrt (x*4)
{-
data Point = Point { _x :: Double, _y :: Double }deriving (Show)
data Atom = Atom { _element :: String, _point :: Point } deriving (Show)
data Molecule = Molecule { _atoms :: [Atom] } deriving (Show)
Control.Lens.makeLenses ''Atom
Control.Lens.makeLenses ''Point
Control.Lens.makeLenses ''Molecule
shiftAtomX :: Atom -> Atom
shiftAtomX = over (point . x) (+ 1)
shiftMoleculeX :: Molecule -> Molecule
shiftMoleculeX = over (atoms . traverse . point . x) (+ 1)
atom1 = Atom { _element = "C", _point = Point { _x = 1.0, _y = 2.0 } }
atom2 = Atom { _element = "O", _point = Point { _x = 3.0, _y = 4.0 } }
molecule = Molecule { _atoms = [atom1, atom2] }
newMolecule = shiftMoleculeX molecule
t = atom1 ^. (point . x)
-}

max_stars_repo_name: rickshilling/game-theory

max_issues_repo_name: rickshilling/game-theory

max_forks_repo_name: rickshilling/game-theory


text:
module Filter (FilterOptions(..), parseOpts, lowPass, convolve) where

import           Control.Monad.State.Strict
import           Data.Coerce
import           Data.Complex
import qualified Data.Vector.Storable as VS
import qualified Data.Vector.Storable.Mutable as VSM
import           Foreign.C.Types
import           Foreign.ForeignPtr
import           Foreign.Ptr
import           Foreign.Storable
import           IQ
import           Options.Applicative
import           Pipes

data FilterOptions = FilterOptions
  { cutoff :: Int
  , transitionWidth :: Int }

parseOpts :: Parser FilterOptions
parseOpts = FilterOptions
         <$> option auto
             ( long "cutoff-frequency"
            <> short 'c'
            <> value 15000
            <> showDefault
            <> metavar "CUTOFF_FREQ"
            <> help "Cutoff frequency, value in Hz")
         <*> option auto
             ( long "transition-width"
            <> short 't'
            <> value 10000
            <> showDefault
            <> metavar "TRANS_WIDTH"
            <> help "Transition width, value in Hz")

calculateNumTaps :: Int -> Int -> Int
calculateNumTaps sr tw = round $ 70.0 * sr' / (22.0 * tw')
  where sr' = fromIntegral sr
        tw' = fromIntegral tw

hann :: Floating a => Int -> [a]
hann n = map (\x -> sin (pi * (fromIntegral x / n')) ** 2) vec
  where vec = [0..n-1]
        n' =  fromIntegral n

-- sin (pi * x) / x. The cutoff frequency controls the sampling frequency of ths
-- sinc function: Sf = 2 * Fc, where Fc is expressed as a fraction of the
-- original signal sampling rate.
sinc :: Floating a =>  Int -> Int -> Int -> [a]
sinc len sr fc = map ((fac *) . sinc') [-(len `div` 2)..(len `div` 2)]
  where sinc' x | x == 0 = 1
                | otherwise =  let n = fromIntegral x in sin (fac * n * pi) / (fac * n * pi)
        fac = 2 * (fromIntegral fc / fromIntegral sr)

lowPass :: Int -> FilterOptions -> VS.Vector Float
lowPass sr opts = VS.fromList x
  where nTaps = calculateNumTaps sr (transitionWidth opts)
        window = hann nTaps
        x = zipWith (*) window $ sinc nTaps sr (cutoff opts)

convolve :: VS.Vector Float -> Pipe IQ IQ IO ()
convolve kernel = do
  let klen = VS.length kernel
      bufSz = 1024
      mkCInt = CInt . fromIntegral

  inBuf <- lift $ VSM.replicate bufSz 0.0
  outBuf <- lift $ VSM.replicate bufSz 0.0
  stateBuf <- lift $ VSM.replicate klen (0.0 :+ 0.0 :: IQ)
  stateHeadPos <- lift mallocForeignPtr
  lift . withForeignPtr stateHeadPos $ flip poke 0

  forever $ do
    forM_ [0..bufSz-1] $ (\i -> await >>= lift . VSM.write inBuf i)

    lift $ VS.unsafeWith (coerce kernel) $ \kBuf ->
      VSM.unsafeWith (coerce stateBuf) $ \stateBuf ->
        VSM.unsafeWith (coerce inBuf) $ \inBuf ->
          VSM.unsafeWith (coerce outBuf) $ \outBuf ->
            withForeignPtr stateHeadPos $ \headPos ->
              conv_c (mkCInt klen) (mkCInt bufSz) headPos kBuf stateBuf inBuf outBuf

    forM_ [0..bufSz-1] $ lift . VSM.read outBuf >=> yield
{-# INLINE convolve #-}

foreign import ccall "convolve" conv_c :: CInt -> CInt -> Ptr CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
module MachineLearning.SoftmaxClassifierTest
(
  tests
)

where

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

import MachineLearning.DataSets (dataset2)

import qualified Numeric.LinearAlgebra as LA
import qualified MachineLearning as ML
import MachineLearning.Optimization
import MachineLearning.SoftmaxClassifier

(x, y) = ML.splitToXY dataset2

model = MultiClass (Softmax 2)

x1 = ML.addBiasDimension x
onesTheta :: LA.Vector LA.R
onesTheta = LA.konst 1 (2 * LA.cols x1)
zeroTheta :: LA.Vector LA.R
zeroTheta = LA.konst 0 (2 * LA.cols x1)

processX muSigma x = ML.addBiasDimension $ ML.featureNormalization muSigma $ ML.mapFeatures 6 x

muSigma = ML.meanStddev (ML.mapFeatures 6 x)
x2 = processX muSigma x

xPredict = LA.matrix 2 [ -0.5, 0.5
                       , 0.2, -0.2
                       , 1, 1
                       , 1, 0
                       , 0, 0
                       , 0, 1]
xPredict2 = processX muSigma xPredict
yExpected = LA.vector [1, 1, 0, 0, 1, 0]

checkSoftmaxGradient theta eps lambda = minimum . take 5 $ map check [eps, eps+0.001 ..]
  where check e = checkGradient model lambda x1 y theta e

gradientCheckingEps :: Double
gradientCheckingEps = 3e-2

eps = 0.000001

initialTheta = LA.konst 0.001 (2 * LA.cols x2)
(thetaGD, optPathGD) = minimize (GradientDescent 0.0005) model eps 150 (L2 1) x2 y initialTheta
(thetaCGFR, optPathCGFR) = minimize (ConjugateGradientFR 0.05 0.2) model eps 30 (L2 1) x2 y initialTheta
(thetaCGPR, optPathCGPR) = minimize (ConjugateGradientPR 0.05 0.3) model eps 30 (L2 1) x2 y initialTheta

showOptPath optPath = show $  (LA.toColumns optPath) !! 1

tests = [  testGroup "gradient checking" [
              testCase "non-zero theta, non-zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 (L2 2))
              , testCase "zero theta, non-zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 (L2 2))
              , testCase "non-zero theta, zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 (L2 0))
              , testCase "zero theta, zero lambda" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 (L2 0))
              , testCase "non-zero theta, no reg" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient onesTheta 1e-3 RegNone)
              , testCase "zero theta, no reg" $ assertApproxEqual "" gradientCheckingEps 0 (checkSoftmaxGradient zeroTheta 1e-3 RegNone)
              ]

        , testGroup "learn" [
            testCase "Gradient Descent" $ assertVector (showOptPath optPathGD) 0.01 yExpected (hypothesis model xPredict2 thetaGD)
            , testCase "Conjugate Gradient FR" $ assertVector (showOptPath optPathCGFR) 0.01 yExpected (hypothesis model xPredict2 thetaCGFR)
            , testCase "Conjugate Gradient PR" $ assertVector (showOptPath optPathCGPR) 0.01 yExpected (hypothesis model xPredict2 thetaCGPR)
            ]
        ]

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
{-# LANGUAGE DataKinds, RankNTypes, TypeFamilies #-}
{-# OPTIONS_GHC -Wno-incomplete-uni-patterns #-}
module TestSingleGradient (testTrees) where

import Prelude

import qualified Data.Array.ShapedS as OS
import qualified Data.Strict.Vector as Data.Vector
import qualified Data.Vector.Generic as V
import           Numeric.LinearAlgebra (Vector)
import           Test.Tasty
import           Test.Tasty.HUnit hiding (assert)

import HordeAd hiding (sumElementsVectorOfDual)
import HordeAd.Core.DualClass (IsPrimal, dAdd, dScale)

import TestCommon

testTrees :: [TestTree]
testTrees = [ testDReverse0
            , testDReverse1
            , testPrintDf
            , testDForward
            , testDFastForward
            , quickCheckForwardAndBackward
            , readmeTests
            , readmeTestsV
            ]

dReverse0
  :: HasDelta r
  => (DualNumberVariables 'DModeGradient r
      -> DualMonadGradient r (DualNumber 'DModeGradient r))
  -> [r]
  -> ([r], r)
dReverse0 f deltaInput =
  let ((results, _, _, _), value) =
        dReverse 1 f (V.fromList deltaInput, V.empty, V.empty, V.empty)
  in (V.toList results, value)

fX :: DualMonad 'DModeGradient Float m
   => DualNumberVariables 'DModeGradient Float
   -> m (DualNumber 'DModeGradient Float)
fX variables = do
  let x = var0 variables 0
  return x

fX1Y :: DualMonad 'DModeGradient Float m
     => DualNumberVariables 'DModeGradient Float
     -> m (DualNumber 'DModeGradient Float)
fX1Y variables = do
  let x = var0 variables 0
      y = var0 variables 1
  x1 <- x +\ 1
  x1 *\ y

fXXY :: DualMonad 'DModeGradient Float m
     => DualNumberVariables 'DModeGradient Float
     -> m (DualNumber 'DModeGradient Float)
fXXY variables = do
  let x = var0 variables 0
      y = var0 variables 1
  xy <- x *\ y
  x *\ xy

fXYplusZ :: DualMonad 'DModeGradient Float m
         => DualNumberVariables 'DModeGradient Float
         -> m (DualNumber 'DModeGradient Float)
fXYplusZ variables = do
  let x = var0 variables 0
      y = var0 variables 1
      z = var0 variables 2
  xy <- x *\ y
  xy +\ z

fXtoY :: DualMonad 'DModeGradient Float m
      => DualNumberVariables 'DModeGradient Float
      -> m (DualNumber 'DModeGradient Float)
fXtoY variables = do
  let x = var0 variables 0
      y = var0 variables 1
  x **\ y

freluX :: DualMonad 'DModeGradient Float m
       => DualNumberVariables 'DModeGradient Float
       -> m (DualNumber 'DModeGradient Float)
freluX variables = do
  let x = var0 variables 0
  reluAct x

testDReverse0 :: TestTree
testDReverse0 = testGroup "Simple dReverse application tests" $
  map (\(txt, f, v, expected) ->
        testCase txt $ dReverse0 f v @?= expected)
    [ ("fX", fX, [99], ([1.0],99.0))
    , ("fX1Y", fX1Y, [3, 2], ([2.0,4.0],8.0))
    , ("fXXY", fXXY, [3, 2], ([12.0,9.0],18.0))
    , ("fXYplusZ", fXYplusZ, [1, 2, 3], ([2.0,1.0,1.0],5.0))
    , ( "fXtoY", fXtoY, [0.00000000000001, 2]
      , ([2.0e-14,-3.2236188e-27],9.9999994e-29) )
    , ("fXtoY2", fXtoY, [1, 2], ([2.0,0.0],1.0))
    , ("freluX", freluX, [-1], ([0.0],0.0))
    , ("freluX2", freluX, [0], ([0.0],0.0))
    , ("freluX3", freluX, [0.0001], ([1.0],1.0e-4))
    , ("freluX4", freluX, [99], ([1.0],99.0))
    , ("fquad", fquad, [2, 3], ([4.0,6.0],18.0))
    , ("scalarSum", vec_omit_scalarSum_aux, [1, 1, 3], ([1.0,1.0,1.0],5.0))
    ]

vec_omit_scalarSum_aux
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
vec_omit_scalarSum_aux vec = returnLet $ foldlDual' (+) 0 vec

sumElementsV
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
sumElementsV variables = do
  let x = var1 variables 0
  returnLet $ sumElements0 x

altSumElementsV
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
altSumElementsV variables = do
  let x = var1 variables 0
  returnLet $ altSumElements0 x

-- hlint would complain about spurious @id@, so we need to define our own.
id2 :: a -> a
id2 x = x

sinKonst
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonst variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    sin x + (id2 $ id2 $ id2 $ konst1 1 2)

sinKonstOut
  :: ( DualMonad d r m
     , Floating (Out (DualNumber d (Vector r))) )
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonstOut variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    unOut $ sin (Out x) + Out (id2 $ id2 $ id2 $ konst1 1 2)

sinDelayed :: (Floating a, IsPrimal d a) => DualNumber d a -> DualNumber d a
sinDelayed (D u u') = delayD (sin u) (dScale (cos u) u')

plusDelayed :: (Floating a, IsPrimal d a)
            => DualNumber d a -> DualNumber d a -> DualNumber d a
plusDelayed (D u u') (D v v') = delayD (u + v) (dAdd u' v')

sinKonstDelay
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonstDelay variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    sinDelayed x `plusDelayed` (id2 $ id2 $ id2 $ konst1 1 2)

powKonst
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
powKonst variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    x ** (sin x + (id2 $ id2 $ id2 $ konst1 (sumElements0 x) 2))

powKonstOut
  :: ( DualMonad d r m
     , Floating (Out (DualNumber d (Vector r))) )
  => DualNumberVariables d r -> m (DualNumber d r)
powKonstOut variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    x ** unOut (sin (Out x)
                + Out (id2 $ id2 $ id2 $ konst1 (sumElements0 x) 2))

powKonstDelay
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
powKonstDelay variables = do
  let x = var1 variables 0
  return $ sumElements0 $
    x ** (sinDelayed x
          `plusDelayed` (id2 $ id2 $ id2 $ konst1 (sumElements0 x) 2))

sinKonstS
  :: forall d r m. DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonstS variables = do
  let x = varS variables 0
  return $ sumElements0 $ fromS1
    ((sin x + (id2 $ id2 $ id2 $ konstS 1))
       :: DualNumber d (OS.Array '[2] r))

sinKonstOutS
  :: forall r d m. ( DualMonad d r m
                   , Floating (Out (DualNumber d (OS.Array '[2] r))) )
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonstOutS variables = do
  let x = varS variables 0
  return $ sumElements0 $ fromS1
    (unOut (sin (Out x) + Out (id2 $ id2 $ id2 $ konstS 1))
       :: DualNumber d (OS.Array '[2] r))

sinKonstDelayS
  :: forall d r m. DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
sinKonstDelayS variables = do
  let x = varS variables 0
  return $ sumElements0 $ fromS1
    ((sinDelayed x `plusDelayed` (id2 $ id2 $ id2 $ konstS 1))
       :: DualNumber d (OS.Array '[2] r))

dReverse1
  :: (r ~ Float, d ~ 'DModeGradient)
  => (DualNumberVariables d r -> DualMonadGradient r (DualNumber d r))
  -> [[r]]
  -> ([[r]], r)
dReverse1 f deltaInput =
  let ((_, results, _, _), value) =
        dReverse 1 f
          (V.empty, V.fromList (map V.fromList deltaInput), V.empty, V.empty)
  in (map V.toList $ V.toList results, value)

testDReverse1 :: TestTree
testDReverse1 = testGroup "Simple dReverse application to vectors tests" $
  map (\(txt, f, v, expected) ->
        testCase txt $ dReverse1 f v @?= expected)
    [ ("sumElementsV", sumElementsV, [[1, 1, 3]], ([[1.0,1.0,1.0]],5.0))
    , ("altSumElementsV", altSumElementsV, [[1, 1, 3]], ([[1.0,1.0,1.0]],5.0))
    , ( "sinKonst", sinKonst, [[1, 3]]
      , ([[0.5403023,-0.9899925]],2.982591) )
    , ( "sinKonstOut", sinKonstOut, [[1, 3]]
      , ([[0.5403023,-0.9899925]],2.982591) )
    , ( "sinKonstDelay", sinKonstDelay, [[1, 3]]
      , ([[0.5403023,-0.9899925]],2.982591) )
    , ( "powKonst", powKonst, [[1, 3]]
      , ([[108.7523,131.60072]],95.58371) )
    , ( "powKonstOut", powKonstOut, [[1, 3]]
      , ([[108.7523,131.60072]],95.58371) )
    , ( "powKonstDelay", powKonstDelay, [[1, 3]]
      , ([[108.7523,131.60072]],95.58371) )
    ]

testPrintDf :: TestTree
testPrintDf = testGroup "Pretty printing test" $
  map (\(txt, f, v, expected) ->
        testCase txt $ prettyPrintDf False f
          (V.empty, V.fromList (map V.fromList v), V.empty, V.empty)
        @?= expected)
    [ ( "sumElementsV", sumElementsV, [[1 :: Float, 1, 3]]
      , unlines
        [ "let0 DeltaId_0 = SumElements0 (Var1 (DeltaId 0)) 3"
        , "in Var0 (DeltaId 0)" ] )
    , ( "altSumElementsV", altSumElementsV, [[1, 1, 3]]
      , unlines
        [ "let0 DeltaId_0 = Add0"
        , "  (Index0 (Var1 (DeltaId 0)) 2 3)"
        , "  (Add0"
        , "     (Index0 (Var1 (DeltaId 0)) 1 3)"
        , "     (Add0 (Index0 (Var1 (DeltaId 0)) 0 3) Zero0))"
        , "in Var0 (DeltaId 0)" ] )
    , ( "sinKonst", sinKonst, [[1, 3]]
      , unlines
        [ "in SumElements0"
        , "  (Add1"
        , "     (Scale1 [ 0.5403023 , -0.9899925 ] (Var1 (DeltaId 0)))"
        , "     (Konst1 Zero0 2))"
        , "  2" ] )
    , ( "sinKonstOut", sinKonstOut, [[1, 3]]
      , unlines
        [ "in SumElements0"
        , "  (Outline1"
        , "     PlusOut"
        , "     [ [ 0.84147096 , 0.14112 ] , [ 1.0 , 1.0 ] ]"
        , "     [ Outline1 SinOut [ [ 1.0 , 3.0 ] ] [ Var1 (DeltaId 0) ]"
        , "     , Konst1 Zero0 2"
        , "     ])"
        , "  2" ] )
    , ( "powKonst", powKonst, [[1, 3]]
      , unlines
        [ "in SumElements0"
        , "  (Add1"
        , "     (Scale1 [ 4.8414707 , 130.56084 ] (Var1 (DeltaId 0)))"
        , "     (Scale1"
        , "        [ 0.0 , 103.91083 ]"
        , "        (Add1"
        , "           (Scale1 [ 0.5403023 , -0.9899925 ] (Var1 (DeltaId 0)))"
        , "           (Konst1 (SumElements0 (Var1 (DeltaId 0)) 2) 2))))"
        , "  2" ] )
    , ( "powKonstOut", powKonstOut, [[1, 3]]
      , unlines
        [ "in SumElements0"
        , "  (Add1"
        , "     (Scale1 [ 4.8414707 , 130.56084 ] (Var1 (DeltaId 0)))"
        , "     (Scale1"
        , "        [ 0.0 , 103.91083 ]"
        , "        (Outline1"
        , "           PlusOut"
        , "           [ [ 0.84147096 , 0.14112 ] , [ 4.0 , 4.0 ] ]"
        , "           [ Outline1 SinOut [ [ 1.0 , 3.0 ] ] [ Var1 (DeltaId 0) ]"
        , "           , Konst1 (SumElements0 (Var1 (DeltaId 0)) 2) 2"
        , "           ])))"
        , "  2" ] )
    ]

testDForward :: TestTree
testDForward =
 testGroup "Simple dForward application tests" $
  map (\(txt, f, v, expected) ->
        let vp = listsToParameters v
        in testCase txt $ dForward f vp vp @?= expected)
    [ ("fquad", fquad, ([2 :: Double, 3], []), (26.0, 18.0))
    , ( "atanReadmeM", atanReadmeM, ([1.1, 2.2, 3.3], [])
      , (7.662345305800865, 4.9375516951604155) )
    , ( "vatanReadmeM", vatanReadmeM, ([], [1.1, 2.2, 3.3])
      , (7.662345305800865, 4.9375516951604155) )
    ]

testDFastForward :: TestTree
testDFastForward =
 testGroup "Simple dFastForward application tests" $
  map (\(txt, f, v, expected) ->
        let vp = listsToParameters v
        in testCase txt $ dFastForward f vp vp @?= expected)
    [ ("fquad", fquad, ([2 :: Double, 3], []), (26.0, 18.0))
    , ( "atanReadmeM", atanReadmeM, ([1.1, 2.2, 3.3], [])
      , (7.662345305800865, 4.9375516951604155) )
    , ( "vatanReadmeM", vatanReadmeM, ([], [1.1, 2.2, 3.3])
      , (7.662345305800865, 4.9375516951604155) )
    ]

-- The formula for comparing derivative and gradient is due to @awf
-- at https://github.com/Mikolaj/horde-ad/issues/15#issuecomment-1063251319
quickCheckForwardAndBackward :: TestTree
quickCheckForwardAndBackward =
  testGroup "Simple QuickCheck of gradient vs derivative vs perturbation"
    [ quickCheckTest0 "fquad" fquad (\(x, y, _z) -> ([x, y], [], [], []))
    , quickCheckTest0 "atanReadmeM" atanReadmeM
             (\(x, y, z) -> ([x, y, z], [], [], []))
    , quickCheckTest0 "vatanReadmeM" vatanReadmeM
             (\(x, y, z) -> ([], [x, y, z], [], []))
    , quickCheckTest0 "sinKonst" sinKonst  -- powKonst NaNs immediately
             (\(x, _, z) -> ([], [x, z], [], []))
    , quickCheckTest0 "sinKonstOut" sinKonstOut
             (\(x, _, z) -> ([], [x, z], [], []))
    , quickCheckTest0 "sinKonstDelay" sinKonstDelay
             (\(x, _, z) -> ([], [x, z], [], []))
    , quickCheckTest0 "sinKonstS" sinKonstS
             (\(x, _, z) -> ([], [], [], [x, z]))
    , quickCheckTest0 "sinKonstOutS" sinKonstOutS
             (\(x, _, z) -> ([], [], [], [x, z]))
    , quickCheckTest0 "sinKonstDelayS" sinKonstDelayS
             (\(x, _, z) -> ([], [], [], [x, z]))
   ]

-- A function that goes from `R^3` to `R^2`, with a representation
-- of the input and the output tuple that is convenient for interfacing
-- with the library.
atanReadmeOriginal :: RealFloat a => a -> a -> a -> Data.Vector.Vector a
atanReadmeOriginal x y z =
  let w = x * sin y
  in V.fromList [atan2 z w, z * x]

-- Here we instantiate the function to dual numbers
-- and add a glue code that selects the function inputs from
-- a uniform representation of objective function parameters
-- represented as delta-variables (`DualNumberVariables`).
atanReadmeVariables
  :: IsScalar d r
  => DualNumberVariables d r -> Data.Vector.Vector (DualNumber d r)
atanReadmeVariables variables =
  let x : y : z : _ = vars variables
  in atanReadmeOriginal x y z

-- According to the paper, to handle functions with non-scalar results,
-- we dot-product them with dt which, for simplicity, we here set
-- to a record containing only ones. We could also apply the dot-product
-- automatically in the library code (though perhaps we should
-- emit a warning too, in case the user just forgot to apply
-- a loss function and that's the only reason the result is not a scalar?).
-- For now, let's perform the dot product in user code.

-- Here is the function for dot product with ones, which is just the sum
-- of elements of a vector.
sumElementsOfDualNumbers
  :: IsScalar d r
  => Data.Vector.Vector (DualNumber d r) -> DualNumber d r
sumElementsOfDualNumbers = V.foldl' (+) 0

-- Here we apply the function.
atanReadmeScalar
  :: IsScalar d r
  => DualNumberVariables d r -> DualNumber d r
atanReadmeScalar = sumElementsOfDualNumbers . atanReadmeVariables

-- Here we introduce a single delta-let binding (`returnLet`) to ensure
-- that if this code is used in a larger context and repeated,
-- no explosion of delta-expressions can happen.
-- If the code above had any repeated non-variable expressions
-- (e.g., if @w@ appeared twice) the user would need to make it monadic
-- and apply @returnLet@ already there.
atanReadmeM
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
atanReadmeM = returnLet . atanReadmeScalar

-- The underscores and empty vectors are placeholders for the vector,
-- matrix and arbitrary tensor components of the parameters tuple,
-- which we here don't use (above we construct a vector output,
-- but it's a vector of scalar parameters, not a single parameter
-- of rank 1).
atanReadmeDReverse :: HasDelta r
                   => Domain0 r -> (Domain0 r, r)
atanReadmeDReverse ds =
  let ((result, _, _, _), value) =
        dReverse 1 atanReadmeM (ds, V.empty, V.empty, V.empty)
  in (result, value)

readmeTests :: TestTree
readmeTests = testGroup "Simple tests for README"
  [ testCase " Float (1.1, 2.2, 3.3)"
    $ atanReadmeDReverse (V.fromList [1.1 :: Float, 2.2, 3.3])
      @?= (V.fromList [3.0715904, 0.18288425, 1.1761366], 4.937552)
  , testCase " Double (1.1, 2.2, 3.3)"
    $ atanReadmeDReverse (V.fromList [1.1 :: Double, 2.2, 3.3])
      @?= ( V.fromList [ 3.071590389300859
                       , 0.18288422990948425
                       , 1.1761365368997136 ]
          , 4.9375516951604155 )
  ]

-- And here's a version of the example that uses vector parameters
-- (quite wasteful in this case) and transforms intermediate results
-- via a primitive differentiable type of vectors instead of inside
-- vectors of primitive differentiable scalars.

vatanReadmeM
  :: DualMonad d r m
  => DualNumberVariables d r -> m (DualNumber d r)
vatanReadmeM variables = do
  let xyzVector = var1 variables 0
      [x, y, z] = map (index0 xyzVector) [0, 1, 2]
      v = seq1 $ atanReadmeOriginal x y z
  returnLet $ sumElements0 v

vatanReadmeDReverse :: HasDelta r
                    => Domain1 r -> (Domain1 r, r)
vatanReadmeDReverse dsV =
  let ((_, result, _, _), value) =
        dReverse 1 vatanReadmeM (V.empty, dsV, V.empty, V.empty)
  in (result, value)

readmeTestsV :: TestTree
readmeTestsV = testGroup "Simple tests of vector-based code for README"
  [ testCase "V Float (1.1, 2.2, 3.3)"
    $ vatanReadmeDReverse (V.singleton $ V.fromList [1.1 :: Float, 2.2, 3.3])
      @?= ( V.singleton $ V.fromList [3.0715904, 0.18288425, 1.1761366]
          , 4.937552 )
  , testCase "V Double (1.1, 2.2, 3.3)"
    $ vatanReadmeDReverse (V.singleton $ V.fromList [1.1 :: Double, 2.2, 3.3])
      @?= ( V.singleton $ V.fromList [ 3.071590389300859
                                     , 0.18288422990948425
                                     , 1.1761365368997136 ]
          , 4.9375516951604155 )
  ]

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
module AoC2021.Day06 where

import           Data.Semigroup        (stimes)
import           Numeric.LinearAlgebra ((#>), (><))
import qualified Numeric.LinearAlgebra as M
import qualified Parsers               as P

-- Linear (in growth days) time algorithm ----------------------------------------------------------
-- See Leslie Matrices (https://en.wikipedia.org/wiki/Leslie_matrix)

leslieMatrix :: M.Matrix Double -- Use double, to use optimized matrix multiplication of hmatrix
leslieMatrix = (9><9) [ 0, 1, 0, 0, 0, 0, 0, 0, 0 -- Count 1 -> 0
                      , 0, 0, 1, 0, 0, 0, 0, 0, 0 -- Count 2 -> 1
                      , 0, 0, 0, 1, 0, 0, 0, 0, 0 -- Count 3 -> 2
                      , 0, 0, 0, 0, 1, 0, 0, 0, 0 -- Count 4 -> 3
                      , 0, 0, 0, 0, 0, 1, 0, 0, 0 -- Count 5 -> 4
                      , 0, 0, 0, 0, 0, 0, 1, 0, 0 -- Count 6 -> 5,
                      , 1, 0, 0, 0, 0, 0, 0, 1, 0 -- Count 7 -> 6, 0 -> 6 (re-spawn, instead of dying)
                      , 0, 0, 0, 0, 0, 0, 0, 0, 1 -- Count 8 -> 7
                      , 1, 0, 0, 0, 0, 0, 0, 0, 0 -- Count 9 -> 8, 0 -> 8 (spawn child with 8)
                      ]

projectNDays :: Int -> M.Matrix Double -> M.Matrix Double
projectNDays = stimes

day0 :: [Int] -> M.Vector Double
day0 = M.vector . zipWith count [0..8] . repeat
  where count x = fromIntegral . length . filter (==x)

solve' :: Int -> [Int] -> Int
solve' days fishes = round $ M.sumElements dayN
  where dayN = projectNDays days leslieMatrix #> day0 fishes
----------------------------------------------------------------------------------------------------

-- Exponential (in growth days) time algorithm -----------------------------------------------------
-- | Calculate the number of fishes spawned by one fish within d days
f :: Int -> Int
f d = (1+) . sum . map (f . childOffset) $ [0 .. nSpawn - 1]
  where
    childOffset i = d-9-(7*i)
    nSpawn = ceiling $ fromIntegral d / 7

-- | Total number of fishes within d days, given a list of fishes with their initial counter
solve :: Int -> [Int] -> Int
solve d = sum . map (f . (d-))
----------------------------------------------------------------------------------------------------

solve01 :: String -> String
solve01 input = case P.runParser P.parseIntegerList () "" input of
  Left e   -> show e
  Right is -> show (solve' 80 is)

-- | Way to inefficient, maybe solve the sum of f?
solve02 :: String -> String
solve02 input = case P.runParser P.parseIntegerList () "" input of
  Left e   -> show e
  Right is -> show (solve' 256 is)


text:
import Graphics.Rendering.Chart.Easy
import Graphics.Rendering.Chart.Backend.Cairo
import Data.Array.IArray
import Numeric.Transform.Fourier.FFT
import Data.Complex
import GHC.Float

type Signal = Array Integer (Complex Double)

signalSize :: (Integer, Integer)
signalSize = (0, 1000)

rect t = array signalSize [(x, if x < t then 1 else 0) | x <- [start..end]] where
  start = fst signalSize
  end = snd signalSize

makeRenderable :: Signal -> [(Double, Double)]
makeRenderable signal = take half . normalize $ map f $ assocs signal where
  f (i, e) = (fromIntegral i, magnitude e)
  half = length signal `div` 2

normalize :: (Fractional t, Integral a) => [(a, t)] -> [(t, t)]
normalize x = map (\(i, e) -> (fromIntegral i / l, e)) x where
  l = fromIntegral $ length x

main = do
  let trans_x = fft $ rect 20
  toFile def "rect_transform.png" $ do
    layout_title .= "Rect fourier transform"
    layout_x_axis . laxis_generate .= scaledAxis def (0, 0.5)
    setColors [opaque blue]
    plot $ line "transform" $ [makeRenderable trans_x]

max_stars_repo_name: SebastianCallh/transform-charts

max_issues_repo_name: SebastianCallh/transform-charts

max_forks_repo_name: SebastianCallh/transform-charts


text:
module ECC.Estimate (estimate,showEstimate) where

import ECC.Types
import qualified Data.Vector.Unboxed as U
import Statistics.Sample (mean)
import Statistics.Resampling (resample, fromResample, Estimator(..))
import Statistics.Resampling.Bootstrap (bootstrapBCA)
import Statistics.Types (Estimate (..), ConfInt, mkCL, ConfInt(..), confidenceLevel)
import System.Random.MWC (create)
import System.Random.MWC
import Numeric

-- Estimate the lower and upper bounds, given a random generator state,
-- a confidence percentage, and a Bit Errors structure.

-- If there are to many samples (as often happens with good error correcting codes),
-- then we cheat, and combine samples.

estimate :: GenIO -> Double -> MessageLength -> BEs -> IO (Maybe (Estimate ConfInt Double))
estimate g confidence m_len bes
  | sumBEs bes == 0 = return Nothing
  | otherwise =
       -- 1000 is the default for criterion
       do resamples <- resample g [Mean] 1000 sampleU
--          print $ U.length $ fromResample $ head $ resamples
--          print resamples
          return $ Just $ head $ bootstrapBCA (mkCL confidence) sampleU {- [Mean] -} resamples
  where
    sample = map (\ be -> be / fromIntegral m_len)
           $ extractBEs bes
    sampleU = U.fromList sample

-- | Show estimate in an understandable format"
showEstimate :: Estimate ConfInt Double -> String
showEstimate est = showEFloat (Just 2) (estPoint est)
                 $ (" +" ++)
                 $ (if confIntLDX estErr == 0 then ("?" ++)
                    else showsPercent ((confIntUDX estErr) / (estPoint est)))
                 $ (" -" ++)
                 $ (if confIntUDX estErr == 0 then ("?" ++)
                    else showsPercent ((confIntLDX estErr) / (estPoint est)))
                 $ (" " ++)
                 $ showsPercent (confidenceLevel (confIntCL estErr))
                 $ ""
  where
    estErr = estError est

showsPercent f = shows (round (100 * f)) . ("%" ++)

-- (1 - (estLowerBound est) / (estPoint est))))

max_stars_repo_name: ku-fpg/ecc-manifold

max_issues_repo_name: ku-fpg/ecc-manifold

max_forks_repo_name: ku-fpg/ecc-manifold


text:
import Framework
import Data.Complex
import Numeric.LinearAlgebra

-- wave: array in either time or frequency domain
-- type Wave = [Double] -- might change to int or unsigned int, still dont know where i will get data from and how it will be formatted
type Wave = Vector R -- Double Vector

-- sample: the target wave and the components used to try and reconstruct it
data Sample = Sample {target::Wave, components::Matrix R}

type Result = Vector R

-- compsToWave :: Result -> [Wave] -> Wave
-- compsToWave ms cs =
--   do
--     let compsList = zipWith multComp ms cs
--     addComps (replicate (length compsList) 0) compsList

-- multComp :: Double -> Wave -> Wave
-- multComp k w = map ((*) k) w

-- addComps :: Wave -> [Wave] -> Wave
-- addComps base [] = base
-- addComps base (w:ws) = addComps (zipWith (+) base w) ws

compsToWave :: Result -> Matrix R -> Wave
compsToWave weights components = components #> weights

test :: (Num a) => a -> a -> a
test q w = q + w

test2 :: (Num a) => a -> a -> a
test2 q w = q * w

main =
  do
    putStrLn (show (evalAlgs [test, test2] [1,2] [[1,4],[2,5],[3,6]]))

-- result: target::wave, predicted::wave or maybe just multiplicities of the components. but what to do about phase?
-- algorithm: input::instance -> predicted::wave/result
-- errcalc: input:result || predicted::wave, target::wave -> some insightful output, IO String maybe?

max_stars_repo_name: ThorntheLetter/wavestimation

max_issues_repo_name: ThorntheLetter/wavestimation

max_forks_repo_name: ThorntheLetter/wavestimation


text:
module Regression where

import MHMonad
import Distr
import Data.Monoid
import Graphics.Rendering.Chart
import Graphics.Rendering.Chart.Backend.Diagrams
import Graphics.Rendering.Chart.State
import Data.Colour
import Data.Colour.Names
import Data.Default.Class
import Control.Lens
import Statistics.Distribution
import Statistics.Distribution.Normal (normalDistr)
import Statistics.Distribution.Beta (betaDistr)

-- A random linear function
linear :: Meas (Double -> Double)
linear =
  do a <- normal 0 3
     b <- normal 0 3
     let f = \x -> a*x + b
     return f

-- Weight a function to some data points
fit :: Meas (a -> Double) -> [(a,Double)] -> Meas (a -> Double)
fit prior dataset = do f <- prior
                       mapM (\(x,y) ->
                               score $ normalPdf (f x) 0.3 y)
                            dataset
                       return f

dataset :: [(Double, Double)]
dataset = [(0,0.6), (1, 0.7), (2,1.2), (3,3.2), (4,6.8), (5, 8.2), (6,8.4)]

example :: Meas (Double -> Double)
example = fit linear dataset

-- Some graphing routines
plotWeightedPoints :: String -> [(Double,Double,Double)] -> Double -> Double -> Double -> Double -> IO ()
plotWeightedPoints filename dataset xmin xmax ymin ymax =
  let plots = map (\(x,y,w) -> toPlot $ plot_points_style .~ filledCircles 4 ((blend w blue red) `withOpacity` w)
                      $ plot_points_values .~ [(x,y)]
                      $ def) dataset in
  let my_layout = layout_plots .~ plots
                $ layout_x_axis .
                  laxis_generate .~ scaledAxis def (xmin,xmax)
                $ layout_y_axis . laxis_generate .~ scaledAxis def (ymin,ymax)
                $ def in
  let graphic =  toRenderable my_layout in
  do
     putStr ("Generating " ++ filename ++ "...")
     renderableToFile def filename graphic;
     putStrLn (" Done!")
     return ()

plotWeightedLines :: String -> [(Double,Double,Double)] -> IO ()
plotWeightedLines filename points =
  let plots = map (\(x,y,w) -> toPlot $ plot_lines_style . line_color .~ (blend w blue red) `withOpacity` (w)
                      $ plot_lines_values .~ [[ (0.0 :: Double,x) , (10.0 :: Double ,y*10+x) ]]
                      $ def) points in
  let my_dots = plot_points_style .~ filledCircles 4 (opaque black)
              $ plot_points_values .~ dataset
              $ def in
  let my_layout = layout_plots .~  (toPlot my_dots : plots)
                $ layout_x_axis .
                  laxis_generate .~ scaledAxis def (-10,10)
                $ layout_y_axis . laxis_generate .~ scaledAxis def (-10,10)
                $ def in
  let graphic =  toRenderable my_layout in
  do
     putStr ("Generating " ++ filename ++ "...")
     renderableToFile def filename graphic;
     putStrLn (" Done!")
     return ()

maxlist :: [Double] -> Double
maxlist (x:xs) = max x (maxlist xs)

minlist :: [Double] -> Double
minlist (x:xs) = min x (minlist xs)

-- Plot the points drawn from weighted samples
test1 = do
        samples <- weightedsamples example
        let xyw's = map (\(f,w) -> (f 0, f 1 - f 0, w)) $ take 10000 samples
        let m = maximum (map (\(_,_,w) -> w) xyw's)
        let xyws = map (\(x,y,w) -> (x,y, max (w/m) (1/255))) xyw's
        putStr $ show $ m
        putStr $ show $ take 20 xyws
        plotWeightedPoints "weightedpoints.svg" xyws (-6) 6 (-6) 6
        --        plotWeightedLines "weightedtest.svg" xyws
        let xyw's = map (\(f,w) -> (f 0, f 1 - f 0, w)) $ take 100000 samples
        let m = maximum (map (\(_,_,w) -> w) xyw's)
        let xyws = filter (\(x,y,_) -> x > -1 && x < 0 && y > 1 && y < 2) $ map (\(x,y,w) -> (x,y, max (w / m) 0.1)) xyw's
        plotWeightedPoints "weightedpointszoom.svg" xyws (-1) 0 1 2

-- Plot the result of MH
test2 = do
        samples <- mh example
        let xyw's = map (\(f,w) -> (f 0, f 1 - f 0, getProduct w)) $ take 10000 samples
        let m = maximum (map (\(_,_,w) -> w) xyw's)
        let xyws = map (\(x,y,w) -> (x,y, max (w / m) 0.05)) xyw's
        plotWeightedPoints "weightedpointsmh.svg" xyws (-6) 6 (-6) 6
        let xyws = filter (\(x,y,_) -> x > -1 && x < 0 && y > 1 && y < 2) $ map (\(x,y,w) -> (x,y, max (w / m) 0.1)) xyw's
        plotWeightedPoints "weightedpointszoommh.svg" xyws (-1) 0 1 2

every n xs = case drop (n-1) xs of
              (y:ys) -> y : every n ys
              [] -> []

-- More graphing routines
-- epsilon: smallest y axis difference to worry about
-- delta: smallest x axis difference to worry about
interesting_points :: (Double -> Double) -> Double -> Double -> Double -> Double -> [Double] -> [Double]
interesting_points f lower upper epsilon delta acc =
  if abs(upper - lower) < delta then acc
  else
    let mid = (upper - lower) / 2 + lower in
    if abs((f(upper) - f(lower)) / 2 + f(lower) - f(mid)) < epsilon
    then acc
    else interesting_points f lower mid epsilon delta (mid : (interesting_points f mid upper epsilon delta acc))

sample_fun f =
--  [ (x, f x) | x <- [(-0.25),(-0.25+0.1)..6.2]]
  let xs = ((-0.25) : (interesting_points f (-0.25) 6.2 0.3 0.001 [6.2])) in
  map (\x -> (x,f x)) xs

plot_funs :: String -> [(Double,Double)] -> [Double -> Double] -> IO ()
plot_funs filename dataset funs =
  let graphs  = map sample_fun funs                 in
  let my_lines  = plot_lines_style . line_color .~ blue `withOpacity` 0.01
                $ plot_lines_values .~ graphs $ def in
  let my_dots = plot_points_style .~ filledCircles 4 (opaque black)
              $ plot_points_values .~ dataset
              $ def in
  let my_layout = layout_plots .~ [toPlot my_lines , toPlot my_dots]
                $ layout_x_axis .
                  laxis_generate .~ scaledAxis def (0,6)
                $ layout_y_axis . laxis_generate .~ scaledAxis def (-2,10)
                $ def in
  let graphic =  toRenderable my_layout in
  do
     putStr ("Generating " ++ filename ++ "...")
     renderableToFile def filename graphic;
     putStrLn (" Done!")
     return ()

-- plot the results of linear regression
test3 =
  do
    fs' <- mh $ fit linear dataset
    let fs = map (\(f,_)->f) $ take 500 $ every 100 $ drop 1000 $ fs'
    plot_funs "lin-reg.svg" dataset fs

-- plot the linear regression prior
test3b =
  do
    fs' <- mh $ linear
    let fs = map (\(f,_)->f) $ take 500 $ every 100 $ drop 1000 $ fs'
    plot_funs "lin-reg-prior.svg" dataset fs

-- poissonPP lower upper rate returns a random list of points between lower and upper coming from a Poisson process with given rate
poissonPP :: Double -> Double -> Double -> Meas [Double]
poissonPP lower upper rate =
  do gap <- exponential rate
     if (lower + gap > upper)
     then return []
     else do rest <- poissonPP (lower + gap) upper rate
             return ((lower + gap) : rest)

-- Splice takes a random list of points, a random function,
-- and returns a new random function that is a piecewise version of the
-- first one.
splice :: Meas [Double] -> Meas (Double -> Double) -> Meas (Double -> Double)
splice pointProcess randomFun =
  do xs <- pointProcess
     fs <- mapM (\_ -> randomFun) xs
     default_f <- randomFun
     let h :: [(Double, Double -> Double)] -> Double -> Double
         h [] x = default_f x
         h ((a,f):xfs) x | x <= a = f x
         h ((a,f):xfs) x | x > a = h xfs x
     return (h (zip xs fs) )

test4 =
  do
    fs' <- mh $ fit (splice (poissonPP 0 6 0.1) linear) dataset
    let fs = map (\(f,_)->f) $ take 500 $ every 1000 $ drop 100000 $ fs'
    plot_funs "piecewise-lin-reg.svg" dataset fs

max_stars_repo_name: PiotrJander/oplss-staton

max_issues_repo_name: PiotrJander/oplss-staton

max_forks_repo_name: PiotrJander/oplss-staton


text:
module Main where
import Plot
import Graphics.GD
import Data.Complex

imageDim :: Size
imageDim = (600,600)

-- How many iterations before  we declare the point doesn't diverge
maxIter = 255

-- Coordinates in the complex plane to graph
-- First coordinate must be top left
mbWindow :: Window
mbWindow = ( (-1.8, 0.02), (-1.7, -0.08) ) -- smaller ship on horizon

main :: IO ()
main = drawPlot drawShip imageDim mbWindow "burning-ship.png"

drawShip :: Coordinate -> Color
drawShip (x,y) = colorIterations $ ship (x :+ y) (0 :+ 0) 0

-- Count number of iterations for the given point to diverge
-- Non-diverging points return 0
ship :: Complex Double -- Coordinate to test
           -> Complex Double -- Iterating Z value
           -> Int -- Current iteration
           -> Int -- Iterations before diverging
ship c z iter
    | iter > maxIter = 0
    | otherwise = let z' = (abs(realPart z) :+ abs(imagPart z))^2 + c in
                  if magnitude z' > 2
                  then iter
                  else ship c z' (iter+1)

colorIterations :: Int -> Color
colorIterations x
    | x > maxIter = rgb 255 255 255
    | otherwise = rgb x x x

max_stars_repo_name: typeclasses/fractal-sets

max_issues_repo_name: typeclasses/fractal-sets

max_forks_repo_name: typeclasses/fractal-sets