text:
{-# LANGUAGE Arrows #-}

module Control.SF.AuxFunctions (
    Event, edge, quantize,
    presentFFT, fftA,
    toMSF, toRealTimeMSF,
    Control, buffer
) where

import Prelude hiding (init)
import Control.Arrow
import Control.CCA.Types

import Numeric.FFT (fft)
import Data.Complex
import Data.Map (Map)
import qualified Data.Map as Map

-- | Alternative for working with Math.FFT instead of Numeric.FFT
--import qualified Math.FFT as FFT
--import Data.Array.IArray
--import Data.Array.CArray
--myFFT n lst = elems $ (FFT.dft) (listArray (0, n-1) lst)

-- For use with SF Conversions
import Control.Monad.Fix
import Control.SF.SF
import Control.SF.MSF

import Control.Concurrent.MonadIO
import Data.IORef.MonadIO
import qualified Data.Sequence as S
import Data.Foldable (toList)
import Control.DeepSeq

--------------------------------------
-- Generic Types and Functions
--------------------------------------

type Event a = Maybe a

edge :: ArrowInit a => a Bool Bool
edge = proc b -> do
    prev <- init False -< b
    returnA -< prev && not b

-- | Scrutinizes n samples at a time, updating after k new values from a signal function
quantize :: ArrowInit a => Int -> Int -> a b (Event [b])
quantize n k = proc d -> do
    rec (ds,c) <- init ([],0) -< (take n (d:ds), c+1)
    returnA -< if c >= n && c `mod` k == 0 then Just ds else Nothing

--------------------------------------
-- Fast Fourier Transform
--------------------------------------

-- | Converts the vector result of a dft into a map from frequency to magnitude.
--   One common use is:
--      fftA >>> arr (fmap $ presentFFT clockRate)
presentFFT :: Double -> [Double] -> Map Double Double
presentFFT clockRate a = Map.fromList $ map mkAssoc (zip [0..(length a)] a) where
    mkAssoc (i,c) = (freq, c) where
        samplesPerPeriod = fromIntegral (length a)
        freq = fromIntegral i * (clockRate / samplesPerPeriod)

-- | Given a quantization frequency (the number of samples between each
--   successive FFT calculation) and a fundamental period, this will decompose
--   the input signal into its constituent frequencies.
--   NOTE: The fundamental period must be a power of two!
fftA :: ArrowInit a => Int -> Int -> a Double (Event [Double])
fftA qf fp = proc d -> do
    carray <- quantize fp qf -< d :+ 0
    returnA -< fmap (map magnitude . take (fp `div` 2) . fft) carray

--------------------------------------
-- Signal Function Conversions
--------------------------------------

-- | The following two functions are for lifting SFs to MSFs.  The first
--   one is a quick and dirty solution, and the second one appropriately
--   converts a simulated time SF into a real time one.
toMSF :: Monad m => SF a b -> MSF m a b
toMSF (SF sf) = MSF h
    where
      h a = return (b, toMSF sf')
        where (b, sf') = sf a

-- | The clockrate is the simulated rate of the input signal function.
--   The buffer is the number of time steps the given signal function is
--   allowed to get ahead of real time.  Thus, the real amount of time
--   that it can get ahead is the buffer divided by the clockrate seconds.
--   The threadHandler is a where the ThreadId of the forked thread is sent.
--
--   The output signal function takes and returns values in real time.
--   The input must be paired with time, and the return values are the
--   list of bs generated in the given time step and a boolean that is
--   true when time is synced and false when the simulation is running
--   slower than real time.  Note that the returned list will be long
--   if the clockrate is much faster than real time and potentially
--   empty if it's slower.
toRealTimeMSF :: (Monad m, MonadIO m, MonadFix m, NFData b) =>
                 Double -> Int -> (ThreadId -> m ()) -> SF a b
              -> MSF m (a, Double) ([b], Bool)
toRealTimeMSF clockrate buffer threadHandler sf = proc (a, t) -> do
    t' <- init 0 -< t
    rec f <- init 0 -< f'
        let tseg = f + t-t'
            (n,f') = properFraction $ tseg * clockrate
    MSF initFun -< (a, n)
  where
      initFun (a, n) = do
        inp <- newIORef a
        out <- newIORef S.empty
        nvar <- newEmptyMVar
        tid <- liftIO $ forkIO $ worker inp out nvar 1 sf
        threadHandler tid
        h inp out nvar (a, n)
      h inp out nvar (a, n) = do
        writeIORef inp a
        if (n > 0)
          then do
            tryPutMVar nvar n
            b <- atomicModifyIORef out (swap . S.splitAt n)
            return ((toList b, S.length b == n), MSF (h inp out nvar))
          else return (([], True),MSF (h inp out nvar))
      worker inp out nvar n (SF sf) = do
        a <- readIORef inp
        let (b, sf') = sf a
        s <- deepseq b $ atomicModifyIORef out (\s -> (s S.|> b, s))
        n' <- if S.length s >= n+buffer then takeMVar nvar else return n
        worker inp out nvar n' sf'
      swap (a,b) = (b,a)

--------------------------------------
-- A Buffering Arrow (beta)
--------------------------------------

data Control = Play | Pause | Record | Dump | Jump Integer
             | ClearEarlier | ClearLater | ClearAll
-- buffer takes a control signal as well as an input stream of data and
-- returns values appropriate to the control signal:
--
buffer :: ArrowInit a => a (Control, b) [b]
buffer = proc (c,x) -> do
    rec s <- init ([],[]) -< s'
        let (ret, s') = case c of
                Play -> next s
                Pause -> ([], s)
                Record -> ([], cons x s)
                Dump -> (llToList s, s)
                Jump n -> ([], nav n s)
                ClearEarlier -> (fst s, ([], snd s))
                ClearLater -> (snd s, (fst s, []))
                ClearAll -> (llToList s, ([],[]))
    returnA -< ret

-- Helper stuff for buffer
-- =======================
type LensList a = ([a],[a])
cons :: a -> LensList a -> LensList a
cons a (e,l) = (a:e, l)

prev :: LensList a -> ([a], LensList a)
prev l@([],_) = ([], l)
prev (a:e, l) = ([a], (e, a:l))

next :: LensList a -> ([a], LensList a)
next l@(_,[]) = ([], l)
next (e, a:l) = ([a], (a:e, l))

nav :: Integer -> LensList a -> LensList a
nav 0 l = l
nav n l | n > 0 = nav (n-1) (snd $ next l)
nav n l | n < 0 = nav (n+1) (snd $ prev l)

llToList :: LensList a -> [a]
llToList (e, l) = e ++ l

max_stars_repo_name: kianwilcox/haskell-music

max_issues_repo_name: kianwilcox/haskell-music

max_forks_repo_name: kianwilcox/haskell-music


text:
{-|
Module      : kMeans
Description : k-Means
Copyright   : (c) Fabrício Olivetti, 2017
License     : GPL-3
Maintainer  : fabricio.olivetti@gmail.com

k-Means Clustering algorithm
-}

module KMeans where

import Data.List
import Data.Ord
--import qualified Numeric.LinearAlgebra as H

import Vector
import Dados

-- |'parseFile' parses a space separated file
-- to a list of lists of Double
parseFile :: String -> [[Double]]
parseFile file = map parseLine (lines file)
  where
    parseLine l = map toDouble (words l)
    toDouble  w = read w :: Double

firstClusters idxs points = [ points !! i | i <- idxs ]

kmeans it points clusters
  | it == 0               = clusters
  | clusters' == clusters = clusters
  | otherwise             = kmeans (it-1) points clusters'
  where
    clusters' = emStep points clusters

emStep points clusters = maximizationStep $ estimationStep points clusters

estimationStep :: [[Double]] -> [[Double]] -> [[[Double]]]
estimationStep points clusters = map unTuple
                               $ groupByValue
                               $ sortByValue
                               $ zip points
                               $ assign points clusters
  where
    unTuple = map fst

maximizationStep :: [[[Double]]] -> [[Double]]
maximizationStep candidates = map mean candidates
  where
    mean xs = (foldl' (.+.) (head xs) (tail xs)) ./ (fromIntegral $ length xs)

assign :: [[Double]] -> [[Double]] -> [Integer]
assign points clusters = map closestTo points
  where
    closestTo p = argmin $ map (euclid p) clusters
    euclid pi ci = sum $ map (^2) $ zipWith (-) pi ci
    argmin xs = fst $ head $ sortByValue $ zip [0..] xs

ext: hs

max_stars_repo_path: 09 - Agrupamento de Dados/KMeans.hs

max_stars_repo_name: douggribeiro/bigdata

max_issues_repo_path: 09 - Agrupamento de Dados/KMeans.hs

max_issues_repo_name: douggribeiro/bigdata

max_forks_repo_path: 09 - Agrupamento de Dados/KMeans.hs

max_forks_repo_name: douggribeiro/bigdata


text:
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE DeriveAnyClass        #-}
{-# LANGUAGE DeriveGeneric         #-}
{-# LANGUAGE FlexibleInstances     #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
{-|
Module      : Grenade.Core.Softmax
Description : Softmax loss layer
Copyright   : (c) Huw Campbell, 2016-2017
License     : BSD2
Stability   : experimental
-}
module Grenade.Layers.Softmax (
    Softmax (..)
  , softmax
  , softmax'
  ) where

import           Data.Serialize

import           Control.DeepSeq              (NFData (..))
import           GHC.Generics                 (Generic)
import           GHC.TypeLits
import           Grenade.Core

import           Numeric.LinearAlgebra.Static as LAS

-- | A Softmax layer
--
--   This layer is like a logit layer, but normalises
--   a set of matricies to be probabilities.
--
--   One can use this layer as the last layer in a network
--   if they need normalised probabilities.
data Softmax = Softmax
  deriving (Show, Generic, NFData)

instance UpdateLayer Softmax where
  type Gradient Softmax = ()
  runUpdate _ _ _ = Softmax

instance RandomLayer Softmax where
  createRandomWith _ _ = return Softmax

instance ( KnownNat i ) => Layer Softmax ('D1 i) ('D1 i) where
  type Tape Softmax ('D1 i) ('D1 i) = S ('D1 i)

  runForwards _ (S1D y) = (S1D y, S1D (softmax y))
  runBackwards _ (S1D y) (S1D dEdy) = ((), S1D (softmax' y dEdy))

instance Serialize Softmax where
  put _ = return ()
  get = return Softmax

softmax :: KnownNat i => LAS.R i -> LAS.R i
softmax xs =
  let xs' = LAS.dvmap exp xs
      s   = LAS.dot xs' 1
  in  LAS.dvmap (/ s) xs'

softmax' :: KnownNat i => LAS.R i -> LAS.R i -> LAS.R i
softmax' x grad =
  let yTy = outer sm sm
      d   = diag sm
      g   = d - yTy
  in  g #> grad
    where
  sm = softmax x

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

instance GNum Softmax where
  _ |* Softmax = Softmax
  _ |+ Softmax = Softmax
  gFromRational _ = Softmax

max_stars_repo_name: koenigmaximilian/grenade

max_issues_repo_name: koenigmaximilian/grenade

max_forks_repo_name: koenigmaximilian/grenade


text:
module Statistics.BBVI.Scheduler
  ( stepTogether
  , stepSeparate
  , unsafeContent
  )
where

import           Control.Monad.ST               ( ST )
import           Data.Propagator
import           Data.Maybe                     ( fromMaybe )
import qualified Data.Vector                   as V
import           Statistics.BBVI.Class
import           Statistics.BBVI.Propagator     ( Time
                                                , time
                                                , DistCell(..)
                                                , DistCellss
                                                , mergeGeneric
                                                , mergeGenericss
                                                )

-- | helper (partial) function to pull content from cell
unsafeContent :: Cell s a -> ST s a
unsafeContent =
  (fromMaybe (error "called unsafe content but no content in cell") <$>)
    . content

-- | update all propagators in one step until first argument reaches
-- quiesence; synonymous with traditional synchronous algorithms
stepTogether :: (a -> b -> ST s (a, b)) -> Cell s a -> Cell s b -> ST s ()
stepTogether f x y = watch x $ \x' -> with y $ \y' -> do
  (upX, upY) <- f x' y'
  -- careful: order below can create a bug. safer would be to use
  -- unsafe content instead of with y
  write y upY
  write x upX

-- | partially optimize each cell with the other held fixed, and
-- repeat until quiesence. Each local optimization will proceed some
-- number of steps, or until the local optimization converges, before
-- triggering the global cell to fire.
stepSeparate
  :: (DistUtil a1, DistUtil a2)
  => Time -- ^ number of local steps
  -> Double
  -> (  DistCell a2
     -> DistCellss a1
     -> ST s (DistCell a2, DistCellss a1)
     ) -- ^ gradient propagator
  -> Cell s (DistCell a2)
  -> Cell s (V.Vector (V.Vector (DistCell a1))) --
  -> ST s ()
stepSeparate nLocal localDelta f x0 ys0 = do
  watch x0 $ \x -> do
    yTemp <- unsafeContent ys0
    ys    <- cellWith $ mergeGenericss
      ((+ nLocal) . time . (V.! 0) . (V.! 0) $ yTemp)
      localDelta
    write ys yTemp
    watch ys $ \ys' -> do
      (_upX, upYs) <- f x ys'
      -- upYs <- f2 x ys'
      write ys upYs
    ysNew <- unsafeContent ys
    write ys0 ysNew
  watch ys0 $ \ys -> do
    xTemp <- unsafeContent x0
    x     <- cellWith $ mergeGeneric (nLocal + time xTemp) localDelta
    write x xTemp
    watch x $ \x' -> do
      (upX, _upYs) <- f x' ys
      -- upX <- f1 x' ys
      write x upX
    xNew <- unsafeContent x
    write x0 xNew

max_stars_repo_name: massma/propagator-bbvi

max_issues_repo_name: massma/propagator-bbvi

max_forks_repo_name: massma/propagator-bbvi


text:
-- ch3 example 3.16

import Statistics.Distribution
import Statistics.Distribution.Hypergeometric

h :: HypergeometricDistribution
h = hypergeometric m l k
  where
    m =  5 --      r
    l = 20 --        N of population size
    k = 10 --          n of samplie size
--                     k elements chosen from
--                   l population, with
--                 m elements of one type and (l-m) of the other

main :: IO ()
main = do
  putStrLn "Among 20, 10 are randomely selected from; what is the probability that the 10 include 5 best."
  putStrLn $ "The syntax is m=r, l=N, k=r, y=k and:"
  putStrLn $ "p 5 = " ++ show (probability h 5)

max_stars_repo_name: raydsameshima/Stat

max_issues_repo_name: raydsameshima/Stat

max_forks_repo_name: raydsameshima/Stat


text:
--
-- Copyright (c) 2019, ERICSSON AB
-- All rights reserved.
--
-- Redistribution and use in source and binary forms, with or without
-- modification, are permitted provided that the following conditions are met:
--
--     * Redistributions of source code must retain the above copyright notice,
--       this list of conditions and the following disclaimer.
--     * Redistributions in binary form must reproduce the above copyright
--       notice, this list of conditions and the following disclaimer in the
--       documentation and/or other materials provided with the distribution.
--     * Neither the name of the ERICSSON AB nor the names of its contributors
--       may be used to endorse or promote products derived from this software
--       without specific prior written permission.
--
-- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-- DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
-- FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
-- DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-- SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-- CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
-- OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
-- OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE StandaloneDeriving #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE ConstraintKinds #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE GADTs #-}
{-# OPTIONS_GHC -Wall -Werror #-}

module Feldspar.Core.Representation
  ( Var(..)
  , VarId(..)
  , TypeF(..)
  , AExpr(..)
  , Info(..)
  , Expr(..)
  , (:->)
  , EqBox(..)
  , Op(..)
  , fvi
  , fviR
  , CBind(..)
  , bvId
  , BindEnv
  , lookupBE
  , extendBE
  , mkLets
  , sharable
  , legalToShare
  , goodToShare
  ) where

import Feldspar.Compiler.Options (Pretty(..))
import Feldspar.Core.NestedTuples (Tuple)
import Feldspar.Core.Types (BoundedInt, Elements, FVal, Length, Index,
                            IntN, IV, MArr, Mut, Par, Size, Type(..),
                            TypeF(..), TypeRep(..))
import Feldspar.Range (Range)

import qualified Data.ByteString.Char8 as B
import Data.Typeable (Typeable, (:~:)(Refl), eqT)
import Data.Hash (Hashable)
import qualified Data.Set as S
import qualified Data.Map.Strict as M

import Data.Bits (Bits, FiniteBits)
import Data.Complex (Complex)
import Data.Ix
import Data.IORef (IORef)

infixr :->
infixl 5 :@
infix 1 :&

-- | Variable identifier
newtype VarId = VarId { varInteger :: Integer }
  deriving (Eq, Ord, Num, Real, Integral, Enum, Ix)

instance Show VarId where
    show (VarId i) = show i

instance Pretty VarId where
  pretty = show

data Var a = Var { varNum :: VarId, varName :: B.ByteString }

instance Eq (Var a) where
  v1 == v2 = varNum v1 == varNum v2

instance Show (Var a) where
  show (Var n _) = "v" ++ show n

-- | Full expressions always have types of the form 'Expr (Full t)' for some 't' to
--   distinguish them from partial applications of operators.
type Full a = a

-- We currently do not use the type constructor :->
type a :-> b = a -> b

-- | The type of information, for instance range information. Currently only size info.
data Info a = Info { infoSize :: Size a }

instance Eq (Size a) => Eq (Info a) where
  Info x == Info y = x == y

instance Show (Size a) => Show (Info a) where
  show (Info x) = show x

-- | Annotated expression, that is, an expression together with extra information,
--   for instance from a program analysis.
data AExpr a where
  (:&) :: TypeF a => {aeInfo :: Info a, aeExpr :: Expr (Full a)} -> AExpr a

deriving instance Eq (AExpr a)

instance Show (AExpr a) where
  show e = showAExpr 0 e ""

showAExpr :: Int -> AExpr a -> String -> String
showAExpr n (i :& (e :: Expr (Full a))) r = '{':inf ++ "} " ++ showExpr n e r
  where inf = show (infoSize i) ++ " : " ++ show (typeRepF :: TypeRep a)

instance Pretty (AExpr a) where
  pretty = show

{- | The main expression type.
     Applications always have an operator at the left end, and are never annotated.

     Note that an operator itself is not a full expression unless its type is of
     the form 'Full t' for some 't'.
-}
data Expr a where
  Sym      :: Typeable a             => Op a -> Expr a
  (:@)     :: Typeable a             => Expr (a -> b) -> AExpr a -> Expr b

instance Show (Expr a) where
  show e = showExpr 0 e ""

showExpr :: Int -> Expr a -> String -> String
showExpr _ (Sym op) r = show op ++ r
showExpr n (f :@ e) r
  = showExpr n f $ "\n" ++ replicate (n + 2) ' ' ++ showAExpr (n + 2) e r

instance Typeable a => Eq (Expr a) where
  Sym op1 == Sym op2 = op1 == op2
  ((f1 :: Expr (a1 -> b1)) :@ e1) == ((f2 :: Expr (a2 -> b2)) :@ e2)
        = case eqT :: Maybe ((a1,b1) :~: (a2,b2)) of
            Nothing -> False
            Just Refl -> f1 == f2 && e1 == e2
  _ == _ = False

-- | A box which makes its contents equal to everything else with the same type
newtype EqBox a = EqBox {unEqBox :: a}

instance Eq (EqBox a) where
  _ == _ = True

instance Show (EqBox a) where
  show _ = "Box"

-- | The main data type for built-in operators as well as let and conditional
--   constructs.
data Op a where
    -- | Array
    Parallel   :: Type a => Op (Length :-> (Index -> a) :-> Full [a])
    Sequential :: (Type a, Type st) =>
                  Op (Length :-> st :-> (Index -> st -> Tuple '[a, st]) :-> Full [a])
    Append     :: Type a => Op ([a] :-> [a] :-> Full [a])
    GetIx      :: Type a => Op ([a] :-> Index :-> Full a)
    SetIx      :: Type a => Op ([a] :-> Index :-> a :-> Full [a])
    GetLength  :: Type a => Op ([a] :-> Full Length)
    SetLength  :: Type a => Op (Length :-> [a] :-> Full [a])

    -- | Binding
    Variable :: Typeable a => Var a -> Op (Full a)
    Lambda  :: Type a => Var a -> Op (b :-> Full (a -> b))
    Let :: Op (a :-> (a -> b) :-> Full b)

    -- | Bits
    BAnd          :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    BOr           :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    BXor          :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    Complement    :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :->       Full a)

    Bit           :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (Index :->       Full a)
    SetBit        :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    ClearBit      :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    ComplementBit :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    TestBit       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full Bool)

    ShiftLU       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    ShiftRU       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    ShiftL        :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> IntN  :-> Full a)
    ShiftR        :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> IntN  :-> Full a)
    RotateLU      :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    RotateRU      :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Index :-> Full a)
    RotateL       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> IntN  :-> Full a)
    RotateR       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> IntN  :-> Full a)
    ReverseBits   :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :->           Full a)

    BitScan       :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Full Index)
    BitCount      :: (Bits a, BoundedInt a, Size a ~ Range a) => Op (a :-> Full Index)

    -- | Complex
    MkComplex :: RealFloat a => Op (a :-> a :-> Full (Complex a))
    RealPart  :: RealFloat a => Op (Complex a :-> Full a)
    ImagPart  :: RealFloat a => Op (Complex a :-> Full a)
    Conjugate :: RealFloat a => Op (Complex a :-> Full (Complex a))
    MkPolar   :: RealFloat a => Op (a :-> a :-> Full (Complex a))
    Magnitude :: RealFloat a => Op (Complex a :-> Full a)
    Phase     :: RealFloat a => Op (Complex a :-> Full a)
    Cis       :: RealFloat a => Op (a :-> Full (Complex a))

    -- | Condition
    Condition  ::                      Op (Bool :-> a :-> a :-> Full a)

    -- | Conversion
    F2I     :: (Integral a, RealFloat b)                     => Op (b :-> Full a)
    -- Type a required for resultType in SizeProp.hs.
    I2N     :: (Type b, Integral a, Num b, Size a ~ Range a) => Op (a :-> Full b)
    B2I     :: (Type a, Integral a)                          => Op (Bool  :-> Full a)
    Round   :: (Integral a, RealFloat b)                     => Op (b :-> Full a)
    Ceiling :: (Integral a, RealFloat b)                     => Op (b :-> Full a)
    Floor   :: (Integral a, RealFloat b)                     => Op (b :-> Full a)

    -- | Elements
    EMaterialize :: Type a => Op (Length :-> Elements a :-> Full [a])
    EWrite       :: Type a => Op (Index :-> a :-> Full (Elements a))
    ESkip        :: Type a => Op (Full (Elements a))
    EPar         :: Type a => Op (Elements a :-> Elements a :-> Full (Elements a))
    EparFor      :: Type a => Op (Length :-> (Index -> Elements a) :-> Full (Elements a))

    -- | Eq
    Equal     :: Eq a => Op (a :-> a :-> Full Bool)
    NotEqual  :: Eq a => Op (a :-> a :-> Full Bool)

    -- | Error
    Undefined :: Type a => Op (Full a)
    Assert    :: Type a => String -> Op (Bool :-> a :-> Full a)

    -- | FFI
    ForeignImport :: String -> EqBox a -> Op (Full a)

    -- | Floating
    Pi      :: Floating a => Op (Full a)
    Exp     :: Floating a => Op (a :-> Full a)
    Sqrt    :: Floating a => Op (a :-> Full a)
    Log     :: Floating a => Op (a :-> Full a)
    Pow     :: Floating a => Op (a :-> a :-> Full a)
    LogBase :: Floating a => Op (a :-> a :-> Full a)
    Sin     :: Floating a => Op (a :-> Full a)
    Tan     :: Floating a => Op (a :-> Full a)
    Cos     :: Floating a => Op (a :-> Full a)
    Asin    :: Floating a => Op (a :-> Full a)
    Atan    :: Floating a => Op (a :-> Full a)
    Acos    :: Floating a => Op (a :-> Full a)
    Sinh    :: Floating a => Op (a :-> Full a)
    Tanh    :: Floating a => Op (a :-> Full a)
    Cosh    :: Floating a => Op (a :-> Full a)
    Asinh   :: Floating a => Op (a :-> Full a)
    Atanh   :: Floating a => Op (a :-> Full a)
    Acosh   :: Floating a => Op (a :-> Full a)

    -- | Fractional
    DivFrac :: Fractional a => Op (a :-> a :-> Full a)

    -- | Future
    MkFuture :: Type a => Op (a :-> Full (FVal a))
    Await    :: Type a => Op (FVal a :-> Full a)

    -- | Integral
    Quot :: (Bits a, Bounded a, Integral a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    Rem  :: (Bits a, Bounded a, Integral a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    Div  :: (Bits a, Bounded a, Integral a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    Mod  :: (Bits a, Bounded a, Integral a, Size a ~ Range a) => Op (a :-> a :-> Full a)
    IExp :: (Bounded a, Integral a, FiniteBits a, Size a ~ Range a) => Op (a :-> a :-> Full a)

    -- | Literal
    Literal :: (Type a, Hashable a) => a -> Op (Full a)

    -- | Logic
    And :: Op (Bool :-> Bool :-> Full Bool)
    Or  :: Op (Bool :-> Bool :-> Full Bool)
    Not :: Op (Bool :->          Full Bool)

    -- | Loop
    ForLoop   :: Type a => Op (Length :-> a :-> (Index -> a -> a) :-> Full a)
    WhileLoop :: Type a => Op (a :-> (a -> Bool) :-> (a -> a) :-> Full a)

    -- | Mutable
    Run :: Type a => Op (Mut a :-> Full a)

    -- | MutableArray
    NewArr    :: Type a => Op (Length :-> a :-> Full (Mut (MArr a)))
    NewArr_   :: Type a => Op (Length :-> Full (Mut (MArr a)))
    GetArr    :: Type a => Op (MArr a :-> Index :-> Full (Mut a))
    SetArr    :: Op (MArr a :-> Index :-> a :-> Full (Mut ()))
    ArrLength :: Op (MArr a :-> Full (Mut Length))

    -- | MutableToPure
    RunMutableArray :: Type a => Op (Mut (MArr a) :-> Full [a])
    WithArray       :: Type b => Op (MArr a :-> ([a] -> Mut b) :-> Full (Mut b))

    -- | MutableReference
    NewRef :: Type a => Op (a :-> Full (Mut (IORef a)))
    GetRef :: Type a => Op (IORef a :-> Full (Mut a))
    SetRef :: Type a => Op (IORef a :-> a :-> Full (Mut ()))
    ModRef :: Type a => Op (IORef a :-> (a -> a) :-> Full (Mut ()))

    -- | Nested tuples
    Cons  :: Type a => Op (a :-> Tuple b :-> Full (Tuple (a ': b)))
    Nil   ::           Op (Full (Tuple '[]))
    Car   :: Type a => Op (Tuple (a ': b) :-> Full a)
    Cdr   ::           Op (Tuple (a ': b) :-> Full (Tuple b))
    Tup   ::           Op (Tuple a :-> Full (Tuple a))

    -- | NoInline
    NoInline :: Type a => Op (a :-> Full a)

    -- | Num
    Abs  :: (Num a, Num (Size a)) => Op (a :-> Full a)
    Sign :: (Num a, Num (Size a)) => Op (a :-> Full a)
    Add  :: (Num a, Num (Size a)) => Op (a :-> a :-> Full a)
    Sub  :: (Num a, Num (Size a)) => Op (a :-> a :-> Full a)
    Mul  :: (Num a, Num (Size a)) => Op (a :-> a :-> Full a)

    -- | Ord
    LTH :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full Bool)
    GTH :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full Bool)
    LTE :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full Bool)
    GTE :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full Bool)
    Min :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full a)
    Max :: (Ord a, Ord (Size a)) => Op (a :-> a :-> Full a)

    -- | Par
    ParRun    :: Type a => Op (Par a :-> Full a)
    ParNew    :: Type a => Op (Full (Par (IV a)))
    ParGet    :: Type a => Op (IV a :-> Full (Par a))
    ParPut    :: Type a => Op (IV a :-> a :-> Full (Par ()))
    ParFork   ::           Op (Par () :-> Full (Par ()))
    ParYield  ::           Op (Full (Par ()))

    -- | RealFloat
    Atan2     :: RealFloat a => Op (a :-> a :-> Full a)

    -- | Save
    Save :: Type a => Op (a :-> Full a)

    -- When are two size prop operators equal?
    -- | SizeProp
    PropSize :: (Type a, Type b) =>
        EqBox (Size a -> Size b) -> Op (a :-> b :-> Full b)

    -- | Switch
    Switch :: Type b => Op (b :-> Full b)

    -- | ConditionM
    ConditionM :: Monad m => Op (Bool :-> m a :-> m a :-> Full (m a))

    -- | LoopM
    While :: (Monad m, Size (m ()) ~ Size ()) => Op (m Bool :-> m a :-> Full (m ()))
    For   :: (Monad m, Size (m ()) ~ Size ()) => Op (Length :-> (Index -> m a) :-> Full (m ()))

    -- | Mutable
    Return :: (Monad m, Size (m a) ~ Size a)         => Op (a    :-> Full (m a))
    Bind   :: (Monad m, Size (m a) ~ Size a)         => Op (m a  :-> (a -> m b) :-> Full (m b))
    Then   :: Monad m                                => Op (m a  :-> m b        :-> Full (m b))
    When   :: Monad m                                => Op (Bool :-> m ()       :-> Full (m ()))

deriving instance Eq (Op a)
deriving instance Show (Op a)

-- | Utility functions

fvi :: AExpr a -> S.Set VarId
fvi (_ :& e) = fviR e

fviR :: Expr a -> S.Set VarId
fviR (Sym (Variable v)) = S.singleton $ varNum v
fviR (Sym (Lambda v) :@ e) = varNum v `S.delete` fvi e
fviR (f :@ e) = fviR f `S.union` fvi e
fviR _ = S.empty

data CBind where
  CBind :: Type a => Var a -> AExpr a -> CBind

instance Eq CBind where
  CBind (v1@Var{} :: Var a) e1 == CBind (v2@Var{} :: Var b) e2
      = case eqT :: Maybe (a :~: b) of
          Nothing -> False
          Just Refl -> v1 == v2 && e1 == e2

instance Show CBind where
  show (CBind v e) = show v ++ " = " ++ show e

bvId :: CBind -> VarId
bvId (CBind v _) = varNum v

mkLets :: ([CBind], AExpr a) -> AExpr a
mkLets (CBind v e1@(Info i1 :& _) : bs, e@(Info i2 :& _))
  = Info i2 :& Sym Let :@ e1 :@ (Info (i1, i2) :& Sym (Lambda v) :@ mkLets (bs, e))
mkLets ([], e) = e

-- | Functions for bind environments
type BindEnv = M.Map VarId CBind

lookupBE :: Typeable a => String -> BindEnv -> Var a -> AExpr a
lookupBE msg bm (v@(Var n _) :: Var a)
               = case M.lookup n bm of
                      Nothing -> error $ msg ++ ": lookupBE does not find variable " ++ show v
                      Just (CBind (Var{} :: Var b) e)
                           -> case eqT :: Maybe (a :~: b) of
                                   Nothing -> error $ msg ++ ": lookupBE finds conflicing types for " ++ show v
                                   Just Refl -> e

extendBE :: BindEnv -> CBind -> BindEnv
extendBE bm b = M.insert (bvId b) b bm

-- | Expressions that can and should be shared
sharable :: AExpr a -> Bool
sharable e = legalToShare e && goodToShare e

-- | Expressions that can be shared without breaking fromCore
legalToShare :: AExpr a -> Bool
legalToShare (_ :& Sym op) = shOp op
legalToShare (_ :& f :@ _)      = shApp f

shApp :: Expr a -> Bool
shApp (f :@ _) = shApp f
shApp (Sym op) = shOp op

shOp :: Op a -> Bool
-- Binding
shOp Variable{} = False
shOp Lambda{}  = False
-- Elements
shOp ESkip     = False
shOp EWrite    = False
shOp EPar      = False
shOp EparFor   = False
-- Monads
shOp Return    = False
shOp Bind      = False
shOp Then      = False
shOp When      = False
-- Monadic arrays
shOp NewArr    = False
shOp NewArr_   = False
shOp GetArr    = False
shOp SetArr    = False
shOp ArrLength = False
-- Monadic loops
shOp For       = False
shOp While     = False
-- MonadRef
shOp NewRef    = False
shOp GetRef    = False
shOp SetRef    = False
shOp ModRef    = False
-- Nested tuples
shOp Cons      = False
shOp Nil       = False
shOp Cdr       = False
-- Everything else
shOp _ = True

-- | Expressions that are expensive enough to be worth sharing
goodToShare :: AExpr a -> Bool
goodToShare (_ :& Sym (Literal (l :: a))) = largeLit (typeRep :: TypeRep a) l
-- The case below avoids constructing a let-binding for an array stored
-- in a tuple. This is beneficial because the select operator is order
-- of magnitudes cheaper than the array copy generated for the let-binding.
-- With a better compilation of array assignments, the need for this
-- special case goes away.
goodToShare (_ :& Sym Car :@ _ :: AExpr a)
  | ArrayType{} <- typeRepF :: TypeRep a
  = False
goodToShare (_ :& _ :@ _) = True
goodToShare _                   = False

largeLit :: TypeRep a -> a -> Bool
largeLit UnitType       _ = False
largeLit BoolType       _ = False
largeLit IntType{}      _ = False
largeLit FloatType      _ = False
largeLit DoubleType     _ = False
largeLit ArrayType{}    l = not $ null l
largeLit ElementsType{} _ = False
largeLit _              _ = True

max_stars_repo_name: Feldspar/feldspar-language

max_issues_repo_name: Feldspar/feldspar-language

max_forks_repo_name: Feldspar/feldspar-language


text:
{-# LANGUAGE RecordWildCards #-}

import Control.Monad.Trans.Except
import Data.Word
import Data.Int
import Data.Complex
import Foreign.C.Types
import Data.Maybe
import Control.Monad
import Data.Monoid

import Control.Error.Util
import Pipes as P
import Pipes.Prelude as P
import Options.Applicative
import Data.Vector.Storable as VS hiding ((++))
import Data.Vector.Generic as VG hiding ((++))
import Graphics.Rendering.OpenGL

import SDR.Util as U
import SDR.RTLSDRStream
import SDR.FFT
import SDR.Plot
import SDR.ArgUtils
import Graphics.DynamicGraph.Waterfall
import Graphics.DynamicGraph.Util

data Options = Options {
    frequency    :: Word32,
    sampleRate   :: Word32,
    gain         :: Maybe Int32,
    fftSize      :: Maybe Int,
    windowWidth  :: Maybe Int,
    windowHeight :: Maybe Int,
    rows         :: Maybe Int,
    colorMap     :: Maybe [GLfloat]
}

parseColorMap :: ReadM [GLfloat]
parseColorMap = eitherReader func
    where
    func "jet"     = return jet
    func "jet_mod" = return jet_mod
    func "hot"     = return hot
    func "bw"      = return bw
    func "wb"      = return wb
    func arg       = Left $ "Cannot parse colour map: `" ++ arg ++ "'"

optParser :: Parser Options
optParser = Options
          <$> option (fmap fromIntegral parseSize) (
                 long "frequency"
              <> short 'f'
              <> metavar "FREQUENCY"
              <> help "Frequency to tune to"
              )
          <*> option (fmap fromIntegral parseSize) (
                 long "samplerate"
              <> short 'r'
              <> metavar "RATE"
              <> help "Sample rate"
              )
          <*> optional (option auto (
                 long "gain"
              <> short 'g'
              <> metavar "GAIN"
              <> help "Tuner gain"
              ))
          <*> optional (option (fmap fromIntegral parseSize) (
                 long "size"
              <> short 's'
              <> metavar "SIZE"
              <> help "FFT bin size. Default is 8192."
              ))
          <*> optional (option auto (
                 long "width"
              <> short 'w'
              <> metavar "WIDTH"
              <> help "Window width. Default is 1024."
              ))
          <*> optional (option auto (
                 long "height"
              <> short 'h'
              <> metavar "HEIGHT"
              <> help "Window height. Default is 480."
              ))
          <*> optional (option auto (
                 long "rows"
              <> short 'r'
              <> metavar "ROWS"
              <> help "Number of rows in waterfall. Default is 1000."
              ))
          <*> optional (option parseColorMap (
                 long "colorMap"
              <> short 'm'
              <> metavar "COLORMAP"
              <> help "Waterfall color map. Default is 'jet_mod'."
              ))

opt :: ParserInfo Options
opt = info (helper <*> optParser) (fullDesc <> progDesc "Draw a dynamic waterall plot of the received spectrum using OpenGL" <> header "RTLSDR Waterfall")

doIt Options{..} = do
    res <- lift setupGLFW
    unless res (throwE "Unable to initilize GLFW")

    let fftSize' =  fromMaybe 8192 fftSize
        window   =  hanning fftSize' :: VS.Vector Double
    str          <- sdrStream ((defaultRTLSDRParams frequency sampleRate) {tunerGain = gain}) 1 (fromIntegral $ fftSize' * 2)
    rfFFT        <- lift $ fftw fftSize'
    rfSpectrum   <- plotWaterfall (fromMaybe 1024 windowWidth) (fromMaybe 480 windowHeight) fftSize' (fromMaybe 1000 rows) (fromMaybe jet_mod colorMap)
    --rfSpectrum   <- plotFill (maybe 1024 id windowWidth) (maybe 480 id windowHeight) fftSize' (maybe jet_mod id colorMap)
    --rfSpectrum   <- plotTexture (maybe 1024 id windowWidth) (maybe 480 id windowHeight) fftSize' fftSize'

    lift $ runEffect $   str
                     >-> P.map (interleavedIQUnsigned256ToFloat :: VS.Vector CUChar -> VS.Vector (Complex Double))
                     >-> P.map (VG.zipWith (flip mult) window . VG.zipWith mult (halfBandUp fftSize'))
                     >-> rfFFT
                     >-> P.map (VG.map ((* (32 / fromIntegral fftSize')) . realToFrac . magnitude))
                     >-> rfSpectrum

main = execParser opt >>= exceptT putStrLn return . doIt

max_stars_repo_name: adamwalker/sdr-apps

max_issues_repo_name: adamwalker/sdr-apps

max_forks_repo_name: adamwalker/sdr-apps


text:
module Formulas where

import           Numeric.LinearAlgebra
import           Control.Monad.IO.Class
import           Control.Monad (replicateM)
import           System.Random
import qualified Data.Vector as V
import           Data.List

type X          = Matrix R
type Beta       = Vector R
type Y          = Vector R
type Block      = V.Vector Int
type Powers     = V.Vector Block
type Input      = V.Vector Double
type Inplist    = V.Vector Input
type Polynomial = (Beta, X)

-- V.Vectors, Vectors and Matrices Transformators
vvToM :: V.Vector (V.Vector Double)
      -> Matrix Double
vvToM = fromLists . V.toList . V.map V.toList

mToVV :: Matrix Double
      -> V.Vector (V.Vector Double)
mToVV = V.fromList . map V.fromList . toLists
-----------------------------------------------------

-- (5) in the article --
f :: Powers
  -> Input
  -> V.Vector Double
f pss inp = V.map oneBlock pss
    where
        oneBlock ps = V.product $ V.zipWith (**) inp $ V.map fromIntegral ps

calcX :: Inplist
      -> Powers -> X
calcX inputs pss = vvToM (V.map (f pss) inputs)

-- (~7) in the article --
cost :: Beta
     -> X
     -> Y
     -> Double
cost b x y = abs $ sumElements $ e**2
    where e = y - x #> b

-- (8) in the article --
normEq :: X
       -> Y
       -> Beta
normEq x y  = (inv (tr x <> x) <> tr x) #> y :: Beta
---------------- Show polynomial ---------------------
showPolynomial :: Beta
               -> Powers
               -> String
showPolynomial coeffs powers =
    intercalate " +" $ zipWith showCoeff coeffs' powers'
        where
            coeffs' = filter isNotZero $ toList coeffs
            powers' = V.toList $ V.map V.toList powers
            showCoeff coeff powerList = " " ++ show_ (round3 coeff) ++ matchVarsPowers powerList
            isNotZero x = round3 x /= 0

matchVarsPowers :: [Int] -> String
-- Example: [1, 4, 2] => x1 * x2^4 * x3^2
matchVarsPowers powers =
    intercalate "" $ zipWith varToPower varNums powers
        where
            varNums = [1..length powers]
            varToPower :: Int -> Int -> String
            varToPower ind 0   = ""
            varToPower ind 1   = " * x" ++ show ind
            varToPower ind pow = " * x" ++ show ind ++ "^" ++ show_ pow

-- show_ :: Num n => n -> n
show_ num = if num > 0 then show num else "(" ++ show num ++ ")"

    -------------------------------------------------------
randomPowerGen :: Int      -- m
               -> Int      -- n
               -> IO (V.Vector (V.Vector Int)) -- m >< n Random Powers
randomPowerGen m n = do
    lists <- replicateM m $ replicateM n $ randomRIO (-4, 4)
    return $ V.fromList $ map V.fromList lists

round3 :: Double -> Double
round3 x = (fromInteger $ round $ x * (10^3)) / (10.0^^3)


text:
{-# LANGUAGE DataKinds, TypeFamilies #-}
module TestConditionalSynth (testTrees) where

import Prelude

import           Data.List (nub, sort, unfoldr)
import qualified Data.Strict.Vector as Data.Vector
import qualified Data.Vector.Generic as V
import qualified Data.Vector.Storable
import           Foreign.Storable (Storable)
import           Foreign.Storable.Tuple ()
import           GHC.Exts (inline)
import           Numeric.LinearAlgebra (Vector)
import qualified Numeric.LinearAlgebra as HM
import           System.Random
import           Test.Tasty
import           Test.Tasty.HUnit hiding (assert)

import HordeAd
import HordeAd.Tool.MnistTools

testTrees :: [TestTree]
testTrees = [ conditionalSynthTests
            ]

-- * A neural net that learns to represent a list of input and result pairs
-- as a function that is a sum of conditionals with linear conditions
-- and linear or zero results.
--
-- The samples are easy: inputs and results are integers and inputs
-- are sorted and not repeating (with input repetition it would not be
-- a function or the output would need to be repeated as well).

bloat :: Int
bloat = 1

lenSynthV :: Int -> Int -> Data.Vector.Vector Int
lenSynthV width nSamples =
  V.fromList $ replicate width (nSamples * 2) ++ [width]
               ++ replicate (bloat * nSamples * 3) width
               ++ replicate 3 (bloat * nSamples)

-- If ai ranges over ps1, bi over ps2 and ci over ps4,
-- the value of the function on input x is
-- the sum of if x * ai + bi > 0 then (x * ai + bi) * ci else 0, which is
-- the sum of if x * ai + bi > 0 then x * ai * ci + bi * ci else 0
-- so each condition depends in a linear way on x and each result,
-- if not 0, depends in a linear way on input x.
--
-- To approximate the samples (a list of input and result pairs on which
-- parameters are trained or tested) using this code, divide the input
-- and multiply result appropriately, see @synthLossSquared@.
synthValue :: forall d r m. DualMonad d r m
           => (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
           -> r
           -> DualNumber d (Vector r)
           -> DualNumber d (Vector r)
           -> DualNumber d (Vector r)
           -> m (DualNumber d r)
synthValue factivation x ps1@(D u _) ps2 ps3 = do
  activated <- factivation $ scale (HM.konst x (V.length u)) ps1 + ps2
  returnLet $ activated <.>! ps3

synthLossSquared :: DualMonad d r m
                 => (DualNumber d (Vector r)
                     -> m (DualNumber d (Vector r)))
                 -> r
                 -> DualNumber d (Vector r)
                 -> DualNumber d (Vector r)
                 -> DualNumber d (Vector r)
                 -> r
                 -> m (DualNumber d r)
synthLossSquared factivation x ps1 ps2 ps3 targ = do
  y <- inline synthValue factivation (x / 1000) ps1 ps2 ps3
  lossSquared (targ / 10000) y  -- smaller target to overcome @tanh@ clamping

-- Inlined to avoid the tiny overhead of calling an unknown function.
sumResultsDual :: forall d r m a. (DualMonad d r m, Storable a)
               => (a -> m (DualNumber d r))
               -> Vector a
               -> m (DualNumber d r)
{-# INLINE sumResultsDual #-}
sumResultsDual f as = do
  let g :: DualNumber d r -> a -> m (DualNumber d r)
      g !acc a = do
        u <- f a
        return $! acc + u
  sumUs <- V.foldM g 0 as
  returnLet sumUs

synthLossAll
  :: forall d r m. DualMonad d r m
  => (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
  -> Data.Vector.Storable.Vector (r, r)
  -> DualNumber d (Vector r)
  -> DualNumber d (Vector r)
  -> DualNumber d (Vector r)
  -> m (DualNumber d r)
synthLossAll factivation samples ps1 ps2 ps3 = do
  let f :: (r, r) -> m (DualNumber d r)
      f (x, y) = inline synthLossSquared factivation x ps1 ps2 ps3 y
  sumResultsDual f samples

sumTrainableInputsS :: forall d r. IsScalar d r
                    => DualNumber d (Vector r)
                    -> Int
                    -> DualNumberVariables d r
                    -> Int
                    -> Data.Vector.Vector (DualNumber d r)
sumTrainableInputsS x offset variables width =
  let f :: Int -> DualNumber d r
      f i = sumTrainableInputsV x (offset + i) variables
  in V.generate width f

splitLayerV :: forall d r m. DualMonad d r m
            => (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
            -> DualNumber d (Vector r)
            -> Int
            -> DualNumberVariables d r
            -> Int
            -> m ( DualNumber d (Vector r)
                 , DualNumber d (Vector r)
                 , DualNumber d (Vector r) )
splitLayerV factivation hiddenVec offset variables width = do
  let multiplied = sumTrainableInputsS hiddenVec offset variables width
      chunkWidth = width `div` 3
      activate :: Int -> m (DualNumber d (Vector r))
      activate n = do
        let v = V.slice (n * chunkWidth) chunkWidth multiplied
        factivation $ seq1 v + var1 variables (offset + width + n)
  a0 <- activate 0
  a1 <- activate 1
  a2 <- activate 2
  return (a0, a1, a2)

synthLossBareTotal
  :: forall d r m. DualMonad d r m
  => (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
  -> (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
  -> (DualNumber d (Vector r) -> m (DualNumber d (Vector r)))
  -> Int
  -> Data.Vector.Storable.Vector (r, r)
  -> DualNumberVariables d r
  -> m (DualNumber d r)
synthLossBareTotal factivation factivationHidden factivationMiddle
                   width samples variables = do
  let (inputs, outputs) = V.unzip samples
      nSamples = V.length samples
      sampleData = inputs <> outputs
      hiddenLayer1 = sumConstantDataL sampleData 0 variables width
                     + var1 variables width  -- bias
  nonlinearLayer1 <- factivationHidden hiddenLayer1
  let offsetMiddle = width + 1
  (ps1, ps2, ps3) <-
    inline splitLayerV factivationMiddle nonlinearLayer1
                       offsetMiddle variables (bloat * nSamples * 3)
  inline synthLossAll factivation samples ps1 ps2 ps3

-- * Tests and generation of random data

-- Pair samples sorted and made unique wrt first element of the pair.
integerPairSamples :: (Storable r, Num r)
                   => (Int, Int) -> Int -> Int
                   -> Data.Vector.Storable.Vector (r, r)
integerPairSamples range seed k =
  let rolls :: RandomGen g => g -> [Int]
      rolls = unfoldr (Just . uniformR range)
      (g1, g2) = split $ mkStdGen seed
      nubInputs :: [Int] -> [Int] -> [Int]
      nubInputs candidates rest =
        let len = length candidates
        in if len == k
           then candidates
           else let (candidatesExtra, restExtra) = splitAt (k - len) rest
                    candidatesUniq = nub $ sort $ candidates ++ candidatesExtra
                in nubInputs candidatesUniq restExtra
      inputs = nubInputs [] (rolls g1)
  in V.zip (V.fromListN k $ map fromIntegral inputs)
           (V.fromListN k $ map fromIntegral $ rolls g2)

gradSmartTestCase
  :: forall r. (HasDelta r, r ~ Double)
  => String
  -> (Int
      -> Data.Vector.Storable.Vector (r, r)
      -> DualNumberVariables 'DModeGradient r
      -> DualMonadGradient r (DualNumber 'DModeGradient r))
  -> Int -> Int -> Int -> Int -> r
  -> TestTree
gradSmartTestCase prefix lossFunction seedSamples
                  nSamples width nIterations expected =
  let makeSamples s =
        integerPairSamples (-1000, 1000) (seedSamples + s) nSamples
      samples = map makeSamples [42, 49 .. 7 * nIterations]
      testSamples = map makeSamples [-1, -2 .. - 100]
      nParams1 = lenSynthV width nSamples
      -- Values from -0.5 to 0.5. TODO: start biases at 1
      params1Init =
        V.imap (\i nPV -> HM.randomVector (33 + nPV + i) HM.Uniform nPV
                          - HM.scalar 0.5)
               nParams1
      parametersInit = (V.empty, params1Init, V.empty, V.empty)
      name = prefix ++ " "
             ++ unwords [ show seedSamples, show nSamples, show width
                        , show (V.length nParams1), show (V.sum nParams1) ]
      f = lossFunction width
  in testCase name $ do
       let (parametersResult, _) =
             sgdAdam f samples parametersInit
                     (initialStateAdam parametersInit)
           (_, values) =
             unzip $ map (\t -> dReverse 1 (f t) parametersResult) testSamples
       (sum values / 100) @?= expected

conditionalSynthTests:: TestTree
conditionalSynthTests = do
 let f = inline synthLossBareTotal reluAct tanhAct tanhAct
 testGroup "synthesizing a sum of linear conditionals matching samples"
  [ gradSmartTestCase "reluAct"
      f 42 10 10  100
      4.740275311294229
  , gradSmartTestCase "reluAct"
      f 42 10 10  10000
      3.83451707827233e-2
  , gradSmartTestCase "reluAct"
      f 42 10 10  100000
      3.135485708489271e-2
  , gradSmartTestCase "reluAct"
      f 42 10 100 100000
      3.2872191198993095e-2
  ]

max_stars_repo_name: Mikolaj/horde-ad

max_issues_repo_name: Mikolaj/horde-ad

max_forks_repo_name: Mikolaj/horde-ad


text:
-- |
-- Module    : Numeric.ApproxEq
-- Copyright : (c) 2011 Aleksey Khudyakov, Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-- Stability   : experimental
-- Portability : portable
--
-- Different implementations of approximate equality for floating
-- point values. There are multiple ways to implement approximate
-- equality. They have different semantics and it's up to programmer
-- to choose right one.
module Numeric.ApproxEq (
    eqRelative
  , eqRelCompl
  , eqAbsolute
  , within
  ) where

import Control.Monad.ST         (runST)

import Data.Complex             (Complex(..))
import Data.Primitive.ByteArray (newByteArray, readByteArray, writeByteArray)
import Data.Int                 (Int64)

-- | Relative difference between two numbers are less than predefined
--   value.  For example 1 is approximately equal to 1.0001 with 1e-4
--   precision. Same is true for 10000 and 10001.
--
--   This method of camparison doesn't work for numbers which are
--   approximately 0. 'eqAbsolute' should be used instead.
eqRelative :: (Fractional a, Ord a)
           => a            -- ^ Relative precision
           -> a -> a -> Bool
eqRelative eps a b = abs (a - b) <= abs (eps * a)
{-# INLINE eqRelative #-}

-- | Relative equality for comlex numbers.
eqRelCompl :: (RealFloat a, Ord a)
           => a                 -- ^ Relative precision
           -> Complex a -> Complex a -> Bool
eqRelCompl eps a b = d <= eps * r
  where
    (d :+ _) = abs (a - b)
    (r :+ _) = abs a
{-# INLINE eqRelCompl #-}

-- | Difference between values is less than specified precision.
eqAbsolute :: (Num a, Ord a)
           => a                 -- ^ Absolute precision
           -> a -> a -> Bool
eqAbsolute eps a b = abs (a - b) <= eps
{-# INLINE eqAbsolute #-}

-- | Compare two 'Double' values for approximate equality, using
-- Dawson's method.
--
-- The required accuracy is specified in ULPs (units of least
-- precision).  If the two numbers differ by the given number of ULPs
-- or less, this function returns @True@.
--
-- Algorithm is based on Bruce Dawson's \"Comparing floating point
-- numbers\":
-- <http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm>
within :: Int                   -- ^ Number of ULPs of accuracy desired.
       -> Double -> Double -> Bool
within ulps a b = runST $ do
  buf <- newByteArray 8
  ai0 <- writeByteArray buf 0 a >> readByteArray buf 0
  bi0 <- writeByteArray buf 0 b >> readByteArray buf 0
  let big  = 0x8000000000000000 :: Int64
      ai | ai0 < 0   = big - ai0
         | otherwise = ai0
      bi | bi0 < 0   = big - bi0
         | otherwise = bi0
  return $! abs (ai - bi) <= fromIntegral ulps

max_stars_repo_name: Shimuuar/numeric-tools

max_issues_repo_name: Shimuuar/numeric-tools

max_forks_repo_name: Shimuuar/numeric-tools


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

module Main ( main ) where

import           Control.Applicative ((<$>))
import qualified Data.Array.Repa     as R
import qualified Data.ByteString     as BS
import           Data.Complex        (Complex (..))
import           Foreign.C.Types     (CDouble, CLLong)
import           Language.J
import           Test.Tasty
import           Test.Tasty.HUnit

main :: IO ()
main = do
#ifdef linux_HOST_OS
    jenv <- jinit libLinux
#else
#ifdef darwin_HOST_OS
    jenv <- jinit (libMac [8,0,7])
#else
#ifdef mingw32_HOST_OS
    jenv <- jinit (libWindows [9,0,1])
#endif
#endif
#endif
    defaultMain $
        testGroup "J dl"
            [ testCase "Performs calculation and has sensible output" (jComp jenv)
            , testCase "Reads back type in the environment" (jType jenv)
            , testCase "Reads a string" (jStr jenv)
            , testCase "Reads a complex array" (jGetC jenv)
            , testCase "Sends an array to J" (jSetA jenv)
            , testCase "Sends a complex array to J" (jSetC jenv)
            , testCase "Uses J to perform a complex calculation" (regress jenv)
            , testCase "Writes strings to J values" (stringRoundtrip jenv)
            -- , testCase "Uses J for something Haskell would have a hard time with" (fill jenv)
            , testCase "Loads a library" (loadNoFail jenv)
            , testCase "Sends an integer array" (printRes jenv)
            ]

loadNoFail :: JEnv -> Assertion
loadNoFail jenv = do
#ifdef linux_HOST_OS
    jLoad jenv (linuxProfile "9.01")
#else
#ifdef darwin_HOST_OS
    jLoad jenv (macProfile [8,0,7])
#else
#ifdef mingw32_HOST_OS
    jLoad jenv (windowsProfile [9,0,1])
#endif
#endif
#endif
    bsDispatch jenv "load'tables/csv'"
    res <- bsOut jenv
    assertBool "Doesn't fail" $
        not ("error: " `BS.isInfixOf` res)

printRes :: JEnv -> Assertion
printRes jenv = do
    setJData jenv "plot_data" $ JIntArr $ R.copyS $ R.map (fromIntegral :: Int -> CLLong) $ R.fromListUnboxed (R.ix1 3) [1,10,2]
    bsDispatch jenv "plot_data"
    res <- bsOut jenv
    res @?= "1 10 2\n"

fill :: JEnv -> Assertion
fill jenv = do
    bsDispatch jenv "random_res =: ? 50 50 $ 1e10"
    res <- getJData jenv "random_res"
    extrExtent res @?= [50, 50]
    where extrExtent :: JData R.DIM2 -> [Int]
          extrExtent (JIntArr res) = R.listOfShape $ R.extent res

regress :: JEnv -> Assertion
regress jenv = do
    let hsArr0 = R.fromListUnboxed (R.ix1 3) [1.0,2.0,3.0]
        hsArr1 = R.fromListUnboxed (R.ix1 3) [2.0,4.0,6.0]
    setJData jenv "xs" (JDoubleArr $ R.copyS $ R.map (realToFrac :: Double -> CDouble) hsArr0)
    setJData jenv "ys" (JDoubleArr $ R.copyS $ R.map (realToFrac :: Double -> CDouble) hsArr1)
    bsDispatch jenv "reg_result =: ys %. xs ^/ i.2"
    res <- getJData jenv "reg_result"
    doubleVect res @?= [5.995204332975845e-15,1.9999999999999971]

stringRoundtrip :: JEnv -> Assertion
stringRoundtrip jenv = do
    setJData jenv "stringy_string" (JString "hello" :: JData R.Z)
    res <- unwrapStr <$> getJData jenv "stringy_string"
    res @?= "hello"

jSetA :: JEnv -> Assertion
jSetA jenv = do
    let hsArr = R.fromListUnboxed (R.ix1 3) [1,3,6]
    setJData jenv "b" (JIntArr $ R.copyS $ R.map (fromIntegral :: Int -> CLLong) hsArr)
    res <- getJData jenv "b"
    intList res @?= [1,3,6]

jSetC :: JEnv -> Assertion
jSetC jenv = do
    let hsArr = R.fromListUnboxed (R.ix1 2) [0.0:+0.0,1.0:+1.0 :: Complex Double]
    let conv = fmap realToFrac :: Complex Double -> Complex CDouble
    setJData jenv "c" (JComplexArr $ R.copyS $ R.map conv hsArr)
    res <- getJData jenv "c"
    complexVect res @?= [0.0:+0.0,1.0:+1.0]

jStr :: JEnv -> Assertion
jStr jenv = do
    bsDispatch jenv "str =: 'hello'"
    res <- getJData jenv "str"
    unwrapStr res @?= "hello"

jGetC :: JEnv -> Assertion
jGetC jenv = do
    bsDispatch jenv "c =: j./~ i.2"
    res <- getJData jenv "c"
    complexVect2 res @?= [0.0:+0.0, 0.0:+1.0, 1.0:+0.0, 1.0:+1.0]

jComp :: JEnv -> Assertion
jComp jenv = do
    bsDispatch jenv "harmonic =: (+/ % #) &.: %"
    bsDispatch jenv "c =: harmonic 1 3 6"
    res <- getJData jenv "c"
    doubleScalar res @?= [2.0]

unwrapStr :: JData R.Z -> BS.ByteString
unwrapStr (JString bs) = bs

doubleVect :: JData R.DIM1 -> [CDouble]
doubleVect (JDoubleArr arr) = R.toList arr

doubleScalar :: JData R.Z -> [CDouble]
doubleScalar (JDoubleArr arr) = R.toList arr

complexVect2 :: JData R.DIM2 -> [Complex CDouble]
complexVect2 (JComplexArr arr) = R.toList arr

complexVect :: JData R.DIM1 -> [Complex CDouble]
complexVect (JComplexArr arr) = R.toList arr

intList :: JData R.DIM1 -> [CLLong]
intList (JIntArr arr) = R.toList arr

jType :: JEnv -> Assertion
jType jenv = do
    bsDispatch jenv "a =: 6?6"
    res <- getJData jenv "a"
    length (intList res) @?= 6

max_stars_repo_name: jitwit/j-hs

max_issues_repo_name: jitwit/j-hs

max_forks_repo_name: jitwit/j-hs


text:
-- An abstract data type for performing multiple linear regression,
-- using HMatrix as the backend for matrix computations.
--
-- types:
--   Matrix a (2D-array)
--   Vector a (vector from the standard vector package)
--
-- Both types are dense, immutable, and strict in all elements and are
-- manipulated as whole blocks.
--
-- matrix product is a <> b, where a and b are matrices
-- matrix-vector product is a #> x, where a is a matrix and x a vector
-- vector-vector dot product is x <.> y
-- the (conjugate) transpose is trans a, where a is a matrix or vector
-- the general linear system solver is <\>
--
-- to construct a 4 by 3 matrix:
-- let a = (4><3) [1,2,3,4,5,6,7,8,9,10,11,12] :: Matrix R
--
-- to construct a vector of double-precision floating points:
-- let x = vector [1,2,3]
--
-- to multiply vectors:
-- x * y
--
-- to add two vectors elementwise:
-- x + y
--
-- to compute vector dot product:
-- x <.> y
--
-- to perform matrix-vector product:
-- m  #> x
--
-- to perform matrix-matrix product
-- m <> n
--
-- to sum elements of a vector or matrix:
-- sumelements m
--

module Regression.HMatrix (
    model_features,
    model_outputs,
    model_weights,
    model_predictions,
    model_rss,
    create_model,
    create_features,
    create_weights,
    predict,
    powers,
    normalize,
    newtons_method,
    newtons_method_with_norm,
    cross_validate,
    gradient_descent
) where

import Numeric.LinearAlgebra hiding (magnitude)
import Numeric.LinearAlgebra.HMatrix hiding (magnitude)
import Data.List hiding (transpose)
import Data.Maybe (fromJust)
import Debug.Trace (trace)
import Stat (range, mean, stdev)
import qualified Data.Vector as V

-- A model is:
-- * a feature matrix with n rows and d columns
-- * a vector of length n representing the observed output
-- * a weight vector containing the calculated vector of optimized weights
-- * the predicted outputs using the optimized weights
-- * the
data Model = MO {
    model_features    :: FeatureMatrix,
    model_outputs     :: FeatureVector,
    model_weights     :: WeightVector,
    model_predictions :: FeatureVector,
    model_rss         :: Double
} deriving (Show)

-- A two-dimensional matrix whose columns contain feature values.
data FeatureMatrix = FM {
    fm_name_indexes:: [(String, Int)],
    fm_values :: Matrix Double
} deriving (Show)

-- A vector of values for a named feature.
data FeatureVector = FV {
    fv_name :: String,
    fv_values :: Vector Double
}
instance Show FeatureVector where
    show (FV name values) =
        name ++ " = " ++ show (toList values)

type Feature a = (String, a -> Double)

type Output a = (String, a -> Double)

type Optimizer = FeatureMatrix -> FeatureVector -> WeightVector

-- A vector of weights.  Each weight corresponds to a feature in the
-- feature matrix used to calculate the weights.
data WeightVector = WV {
    wv_name_indexes :: [(String, Int)],
    wv_values :: Vector Double
}

instance Show WeightVector where
    show (WV name_indexes values) =
        concatMap showWeight (zip name_indexes (toList values))
        where showWeight ((name, i), value) = name ++ " = " ++ show value ++ "\n"

-- Creates a model from a list of records, a list of features, and the output
-- accessor function.  This computes the weight vector as well.  This function
-- does not scale the features nor does it add an intercept feature (a feature
-- whose values are all 1).
create_model :: [a] -> [Feature a] -> Output a -> Optimizer -> Model
create_model rows features (output_name, output) optimizer =
    let fmat         = create_features features rows
        nn           = length rows
        observations = FV output_name (fromList (Data.List.map output rows))
        weights      = optimizer fmat observations
        predictions  = predict fmat weights
        residuals    = rss observations predictions
    in MO fmat observations weights predictions residuals

-- Multiplies a feature matrix by a weights vector to obtain a prediction of
-- output.
predict :: FeatureMatrix -> WeightVector -> FeatureVector
predict (FM _ h) (WV _ w) = FV "" (h #> w)

-- Generates the feature matrix, usually denoted as H, with N rows and
-- D features where D is the length of the feature list.
create_features :: [Feature a] -> [a] -> FeatureMatrix
create_features hs inputs =
    let n = length inputs
        d = length hs
        names = (Data.List.map fst hs)
        name_indexes = Prelude.zip names [0..]
        dat = [h(row) | row <- inputs, (_,h) <- hs]
        h = (n><d) dat
    in FM name_indexes h

-- Computes the residual sum of squares from the observed output and
-- the predicted output.
rss :: (Num a) => FeatureVector -> FeatureVector -> Double
rss (FV _ v1) (FV _ v2) = rss' v1 v2

rss' :: Vector Double -> Vector Double -> Double
rss' v1 v2 = diff <.> diff
    where diff = v1 - v2

-- Takes a list of numbers and turns them into a vector of weights.
create_weights :: FeatureMatrix -> [Double] -> WeightVector
create_weights (FM name_indexes h) weights = WV name_indexes (vector weights)

newtons_method :: Optimizer
newtons_method h y = newtons_method_with_norm 0 h y

newtons_method_with_norm :: Matrix Double -> Optimizer
newtons_method_with_norm l (FM n h) (FV o y) = WV n w
    where w = newtons_method_with_norm' l h y

-- Performs newtons method, finding an estimate of the minimum.  This method
-- uses the formula:
-- w_hat = (H^t * H - lambda*I)^-1 * H^T * y
newtons_method_with_norm' :: Matrix Double -> Matrix Double -> Vector Double -> Vector Double
newtons_method_with_norm' l h y = ((inv (hth - mod)) <> th) #> y
    where th = tr h
          hth = th <> h
          (d,_) = size hth
          mod = l * (ident d) :: Matrix Double

-- Performs gradient descent, updating the weights until the residual
-- sum of squares is less than the epsilon value e, at which point the
-- weight matrix is returned.  n is the step size.
gradient_descent :: Double -> Double -> WeightVector -> Optimizer
gradient_descent e n (WV name w) (FM _ f) (FV _ o) =
    let ft      = tr f
        weights = gradient_descent' e n w f ft o
    in WV name weights

gradient_descent' :: Double -> Double -> Vector Double -> Matrix Double ->
    Matrix Double -> Vector Double -> Vector Double
gradient_descent' e n w h ht y =
    let grad = gradient h ht y w -- (-2H^t(y-Hw))
        grad_len = magnitude grad  -- grad RSS(w) == ||2H^t(y-HW)||
    --in if grad_len < e
    in if (trace ("gradient = " ++ show grad_len) grad_len) < e
        then w
        else let delta = cmap (*(-n)) grad -- (2nH^t(y-Hw))
                 w' = w + delta
             in gradient_descent' e n w' h ht y

-- Compute the magnitude of the given vector.
magnitude :: Vector Double -> Double
magnitude vec = sumElements $ cmap (\y -> y^2) vec

-- Calculates the gradient of the residual sum of squares (-2H^t(y-Hw)).
-- This is used to compute the magnitude of the gradient, to see if the
-- function is minimized.  It is also used to update the weights of the
-- features.
gradient :: Matrix Double -> Matrix Double -> Vector Double ->
    Vector Double -> Vector Double
gradient h ht y w =
    let yhat = h #> w
        err = y - yhat
        prod = ht #> err
        grad = cmap (*(-2)) prod
    in grad

-- Given a feature, computes polynomial powers of that feature from
-- 1 (the original feature) up to and including n.
powers :: Feature a -> Int -> [Feature a]
powers (name, f) n = Data.List.map  (\i -> (name ++ (show i), (\a -> (f a)^i))) [1..n]

-- given a list of values and a feature, return a normalized version of that
-- feature where the mean of the feature is subtracted from each feature and
-- the result is divided by the standard deviation of the feature.
normalize :: [a] -> (a -> Double) -> a -> Double
normalize xs f = (/sdxs) . (+(-meanxs)) . f
    where meanxs = mean (Data.List.map f xs)
          sdxs = stdev (Data.List.map f xs)

-- Compute the average k-fold cross-validation error for a given L2 penalty and
-- a given split size k.
-- Note that the rows should be shuffled before they are provided.
cross_validate :: Int -> [a] -> [Feature a] -> Output a  -> Optimizer -> Double
cross_validate k rows features o@(output_name, output) optimizer = kmean
    where ksplits         = ksplit k rows
          training_data   = map fst ksplits
          valid_data      = map snd ksplits
          trained_models  = map (\rows -> create_model rows features o optimizer) training_data
          trained_weights = map model_weights trained_models
          valid_features  = map (create_features features) valid_data
          predictions     = zipWith (\h w -> predict h w) valid_features trained_weights
          observations    = map (\dat -> FV output_name (fromList (Data.List.map output dat))) valid_data
          valid_rss       = zipWith (\o p -> rss o p) observations predictions
          kmean           = mean valid_rss

-- prepare splits of the data for k-means where k is the number of
-- clusters and the number of values in the test set (the second list
-- in the pair) is n/k where n is the size of the input list.
ksplit :: Int -> [a] -> [([a], [a])]
ksplit k xs = map (ksplit' k xs) [0..k-1]

-- prepare a k split ([a], [a]) for the chunk at index i with list xs
-- and with number of chunks k
ksplit' :: Int -> [a] -> Int -> ([a], [a])
ksplit' k xs i = ((take (sz*i) xs) ++ (drop (sz*(i+1)) xs), take sz (drop (sz*i) xs))
    where sz = length xs `div` k

max_stars_repo_name: markrgrant/regression

max_issues_repo_name: markrgrant/regression

max_forks_repo_name: markrgrant/regression