text:
-- |
-- Module:      Data.Euclidean
-- Copyright:   (c) 2019 Andrew Lelechenko
-- Licence:     BSD3
-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>
--

{-# LANGUAGE CPP                        #-}
{-# LANGUAGE DefaultSignatures          #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MagicHash                  #-}

module Data.Euclidean
  ( Euclidean(..)
  , Field
  , GcdDomain(..)
  , WrappedIntegral(..)
  , WrappedFractional(..)
  , gcdExt
  ) where

import Prelude hiding (quotRem, quot, rem, divMod, div, mod, gcd, lcm, negate, (*), Int, Word)
import qualified Prelude as P
import Control.Exception (throw, ArithException(..))
import Data.Bits (Bits)
import Data.Complex (Complex(..))
import Data.Int (Int, Int8, Int16, Int32, Int64)
import Data.Maybe (isJust)
import Data.Ratio (Ratio)
import Data.Semiring (Semiring(..), Ring(..), (*), minus, isZero, Mod2)
import Data.Word (Word, Word8, Word16, Word32, Word64)
import Foreign.C.Types (CFloat, CDouble)
import GHC.Exts (Int(..), Word(..))
import GHC.Integer.GMP.Internals (gcdInt, gcdWord, gcdInteger, lcmInteger)

import Numeric.Natural

---------------------------------------------------------------------
-- Classes
---------------------------------------------------------------------

-- | 'GcdDomain' represents a
-- <https://en.wikipedia.org/wiki/GCD_domain GCD domain>.
-- This is a domain, where GCD can be defined,
-- but which does not necessarily allow a well-behaved
-- division with remainder (as in 'Euclidean' domains).
--
-- For example, there is no way to define 'rem' over
-- polynomials with integer coefficients such that
-- remainder is always "smaller" than divisor. However,
-- 'gcd' is still definable, just not by means of
-- Euclidean algorithm.
--
-- All methods of 'GcdDomain' have default implementations
-- in terms of 'Euclidean'. So most of the time
-- it is enough to write:
--
-- > instance GcdDomain Foo
-- > instance Euclidean Foo where
-- >   quotRem = ...
-- >   degree  = ...
class Semiring a => GcdDomain a where
  -- | Division without remainder.
  --
  -- prop> \x y -> (x * y) `divide` y == Just x
  -- prop> \x y -> maybe True (\z -> x == z * y) (x `divide` y)
  divide :: a -> a -> Maybe a

  default divide :: (Eq a, Euclidean a) => a -> a -> Maybe a
  divide x y = let (q, r) = quotRem x y in
    if isZero r then Just q else Nothing

  -- | Greatest common divisor. Must satisfy
  --
  -- prop> \x y -> isJust (x `divide` gcd x y) && isJust (y `divide` gcd x y)
  -- prop> \x y z -> isJust (gcd (x * z) (y * z) `divide` z)
  gcd :: a -> a -> a

  default gcd :: (Eq a, Euclidean a) => a -> a -> a
  gcd a b
    | isZero b  = a
    | otherwise = gcd b (a `rem` b)

  -- | Lowest common multiple. Must satisfy
  --
  -- prop> \x y -> isJust (lcm x y `divide` x) && isJust (lcm x y `divide` y)
  -- prop> \x y z -> isNothing (z `divide` x) || isNothing (z `divide` y) || isJust (z `divide` lcm x y)
  lcm :: a -> a -> a

  default lcm :: Eq a => a -> a -> a
  lcm a b
    | isZero a || isZero b = zero
    | otherwise = case a `divide` gcd a b of
      Nothing -> error "lcm: violated gcd invariant"
      Just c  -> c * b

  -- | Test whether two arguments are
  -- <https://en.wikipedia.org/wiki/Coprime_integers coprime>.
  -- Must match its default definition:
  --
  -- prop> \x y -> coprime x y == isJust (1 `divide` gcd x y)
  coprime :: a -> a -> Bool

  default coprime :: Eq a => a -> a -> Bool
  coprime x y = isJust (one `divide` gcd x y)

infixl 7 `divide`

-- | Informally speaking, 'Euclidean' is a superclass of 'Integral',
-- lacking 'toInteger', which allows to define division with remainder
-- for a wider range of types, e. g., complex integers
-- and polynomials with rational coefficients.
--
-- 'Euclidean' represents a
-- <https://en.wikipedia.org/wiki/Euclidean_domain Euclidean domain>
-- endowed by a given Euclidean function 'degree'.
--
-- No particular rounding behaviour is expected of 'quotRem'. E. g.,
-- it is not guaranteed to truncate towards zero or towards negative
-- infinity (cf. 'P.divMod'), and remainders are not guaranteed to be non-negative.
-- For a faithful representation of residue classes one can use
-- <http://hackage.haskell.org/package/mod mod> package instead.
class GcdDomain a => Euclidean a where
  {-# MINIMAL (quotRem | quot, rem), degree #-}
  -- | Division with remainder.
  --
  -- prop> \x y -> y == 0 || let (q, r) = x `quotRem` y in x == q * y + r
  quotRem :: a -> a -> (a, a)
  quotRem x y = (quot x y, rem x y)

  -- | Division. Must match its default definition:
  --
  -- prop> \x y -> quot x y == fst (quotRem x y)
  quot :: a -> a -> a
  quot x y = fst (quotRem x y)

  -- | Remainder. Must match its default definition:
  --
  -- prop> \x y -> rem x y == snd (quotRem x y)
  rem :: a -> a -> a
  rem x y = snd (quotRem x y)

  -- | Euclidean (aka degree, valuation, gauge, norm) function on @a@. Usually @'fromIntegral' '.' 'abs'@.
  --
  -- 'degree' is rarely used by itself. Its purpose
  -- is to provide an evidence of soundness of 'quotRem'
  -- by testing the following property:
  --
  -- prop> \x y -> y == 0 || let (q, r) = x `quotRem` y in (r == 0 || degree r < degree y)
  degree :: a -> Natural

infixl 7 `quot`
infixl 7 `rem`

coprimeIntegral :: Integral a => a -> a -> Bool
coprimeIntegral x y = (odd x || odd y) && P.gcd x y == 1

-- | Execute the extended Euclidean algorithm.
-- For elements @a@ and @b@, compute their greatest common divisor @g@
-- and the coefficient @s@ satisfying @as + bt = g@ for some @t@.
gcdExt :: (Eq a, Euclidean a, Ring a) => a -> a -> (a, a)
gcdExt = go one zero
  where
    go s s' r r'
      | r' == zero = (r, s)
      | otherwise  = case quotRem r r' of
        (q, r'') -> go s' (minus s (times q s')) r' r''
{-# INLINABLE gcdExt #-}

-- | 'Field' represents a
-- <https://en.wikipedia.org/wiki/Field_(mathematics) field>,
-- a ring with a multiplicative inverse for any non-zero element.
class (Euclidean a, Ring a) => Field a

---------------------------------------------------------------------
-- Instances
---------------------------------------------------------------------

instance GcdDomain () where
  divide  = const $ const (Just ())
  gcd     = const $ const ()
  lcm     = const $ const ()
  coprime = const $ const True

instance Euclidean () where
  degree  = const 0
  quotRem = const $ const ((), ())
  quot    = const $ const ()
  rem     = const $ const ()

instance Field ()

instance GcdDomain Mod2 where

instance Euclidean Mod2 where
  degree = const 0
  quotRem x y
    | isZero y  = throw DivideByZero
    | otherwise = (x, zero)

instance Field Mod2

-- | Wrapper around 'Integral' with 'GcdDomain'
-- and 'Euclidean' instances.
newtype WrappedIntegral a = WrapIntegral { unwrapIntegral :: a }
  deriving (Eq, Ord, Show, Num, Integral, Real, Enum, Bits)

instance Num a => Semiring (WrappedIntegral a) where
  plus  = (P.+)
  zero  = 0
  times = (P.*)
  one   = 1
  fromNatural = P.fromIntegral

instance Num a => Ring (WrappedIntegral a) where
  negate = P.negate

instance Integral a => GcdDomain (WrappedIntegral a) where
  divide x y = case x `P.quotRem` y of (q, 0) -> Just q; _ -> Nothing
  gcd     = P.gcd
  lcm     = P.lcm
  coprime = coprimeIntegral

instance Integral a => Euclidean (WrappedIntegral a) where
  degree  = P.fromIntegral . abs . unwrapIntegral
  quotRem = P.quotRem
  quot    = P.quot
  rem     = P.rem

instance GcdDomain Int where
  divide x y = case x `P.quotRem` y of (q, 0) -> Just q; _ -> Nothing
#if MIN_VERSION_integer_gmp(0,5,1)
  gcd (I# x) (I# y) = I# (gcdInt x y)
#else
  gcd     = P.gcd
#endif
  lcm     = P.lcm
  coprime = coprimeIntegral

instance GcdDomain Word where
  divide x y = case x `P.quotRem` y of (q, 0) -> Just q; _ -> Nothing
#if MIN_VERSION_integer_gmp(1,0,0)
  gcd (W# x) (W# y) = W# (gcdWord x y)
#else
  gcd     = P.gcd
#endif
  lcm     = P.lcm
  coprime = coprimeIntegral

instance GcdDomain Integer where
  divide x y = case x `P.quotRem` y of (q, 0) -> Just q; _ -> Nothing
  gcd     = gcdInteger
  lcm     = lcmInteger
  coprime = coprimeIntegral

#define deriveGcdDomain(ty)     \
instance GcdDomain (ty) where { \
;  divide x y = case x `P.quotRem` y of { (q, 0) -> Just q; _ -> Nothing } \
;  gcd     = P.gcd              \
;  lcm     = P.lcm              \
;  coprime = coprimeIntegral    \
}

deriveGcdDomain(Int8)
deriveGcdDomain(Int16)
deriveGcdDomain(Int32)
deriveGcdDomain(Int64)
deriveGcdDomain(Word8)
deriveGcdDomain(Word16)
deriveGcdDomain(Word32)
deriveGcdDomain(Word64)
deriveGcdDomain(Natural)

#define deriveEuclidean(ty)       \
instance Euclidean (ty) where {   \
;  degree  = P.fromIntegral . abs \
;  quotRem = P.quotRem            \
;  quot    = P.quot               \
;  rem     = P.rem                \
}

deriveEuclidean(Int)
deriveEuclidean(Int8)
deriveEuclidean(Int16)
deriveEuclidean(Int32)
deriveEuclidean(Int64)
deriveEuclidean(Integer)
deriveEuclidean(Word)
deriveEuclidean(Word8)
deriveEuclidean(Word16)
deriveEuclidean(Word32)
deriveEuclidean(Word64)
deriveEuclidean(Natural)

-- | Wrapper around 'Fractional'
-- with trivial 'GcdDomain'
-- and 'Euclidean' instances.
newtype WrappedFractional a = WrapFractional { unwrapFractional :: a }
  deriving (Eq, Ord, Show, Num, Fractional)

instance Num a => Semiring (WrappedFractional a) where
  plus  = (P.+)
  zero  = 0
  times = (P.*)
  one   = 1
  fromNatural = P.fromIntegral

instance Num a => Ring (WrappedFractional a) where
  negate = P.negate

instance Fractional a => GcdDomain (WrappedFractional a) where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Fractional a => Euclidean (WrappedFractional a) where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Fractional a => Field (WrappedFractional a)

instance Integral a => GcdDomain (Ratio a) where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Integral a => Euclidean (Ratio a) where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Integral a => Field (Ratio a)

instance GcdDomain Float where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Euclidean Float where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Field Float

instance GcdDomain Double where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Euclidean Double where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Field Double

instance GcdDomain CFloat where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Euclidean CFloat where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Field CFloat

instance GcdDomain CDouble where
  divide x y = Just (x / y)
  gcd        = const $ const 1
  lcm        = const $ const 1
  coprime    = const $ const True

instance Euclidean CDouble where
  degree      = const 0
  quotRem x y = (x / y, 0)
  quot        = (/)
  rem         = const $ const 0

instance Field CDouble

conjQuotAbs :: Field a => Complex a -> Complex a
conjQuotAbs (x :+ y) = x `quot` norm :+ (negate y) `quot` norm
  where
    norm = (x `times` x) `plus` (y `times` y)

instance Field a => GcdDomain (Complex a) where
  divide x y = Just (x `times` conjQuotAbs y)
  gcd        = const $ const one
  lcm        = const $ const one
  coprime    = const $ const True

instance Field a => Euclidean (Complex a) where
  degree      = const 0
  quotRem x y = (quot x y, zero)
  quot x y    = x `times` conjQuotAbs y
  rem         = const $ const zero

instance Field a => Field (Complex a)

max_stars_repo_name: incertia/semirings

max_issues_repo_name: incertia/semirings

max_forks_repo_name: incertia/semirings


text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE Strict           #-}
module STC.Plan
  ( module DFT.Plan
  , module STC.Plan
  ) where

import           Control.Monad        as M
import           Data.Array.Repa      as R
import           Data.Complex
import           Data.List            as L
import           Data.Vector.Storable as VS
import           DFT.Plan
import           System.FilePath
import           System.Random
import           Control.Concurrent.Async

{-# INLINE generateRadomVector #-}
generateRadomVector ::
     (Storable e, Num e, Random e)
  => Int
  -> IO (VS.Vector (Complex e))
generateRadomVector len =
  newStdGen >>= return . VS.map (:+ 0) . VS.fromList . L.take len . randoms

makePlanDiscrete :: FilePath -> DFTPlan -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> IO DFTPlan
makePlanDiscrete folderPath initPlan nx ny numR2Freq numThetaFreq numRFreq numPhiFreq numRhoFreq = do
  importFFTWWisdom (folderPath </> "fftwwisdom.dat")
  let lens =
        [ numThetaFreq * numR2Freq ^ 2
        , nx * ny * numThetaFreq * numRFreq
        , numR2Freq ^ 2
        , numThetaFreq * numRFreq * numPhiFreq * numRhoFreq
        , numRhoFreq * numRFreq
        , numThetaFreq * numR2Freq ^ 2
        , numRFreq
        ]
  (initVec1:initVec2:initVec3:initVec4:initVec5:initVec6:initVec7:_) <-
   mapConcurrently generateRadomVector lens
  lock <- getFFTWLock
  plan <-
    fst <$>
    (dft1dGPlan
       lock
       initPlan
       [numThetaFreq, numR2Freq, numR2Freq]
       [1, 2]
       initVec1 >>= \(plan, vec) ->
       idft1dGPlan
         lock
         plan
         [numThetaFreq, numR2Freq, numR2Freq]
         [1, 2]
         vec >>= \(plan, _) ->
         dft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] initVec2 >>= \(plan, vec) ->
           idft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] vec >>= \(plan, vec) ->
             dft1dGPlan
               lock
               plan
               [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
               [0]
               vec >>= \(plan, vec) ->
               idft1dGPlan
                 lock
                 plan
                 [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
                 [0]
                 vec
            >>= \(plan, _) ->
             dft1dGPlan lock plan [numR2Freq, numR2Freq] [0, 1] initVec3 >>= \(plan, vec) ->
               idft1dGPlan lock plan [numR2Freq, numR2Freq] [0, 1] vec >>= \(plan, _) ->
                 dft1dGPlan
                   lock
                   plan
                   [numRFreq, numRhoFreq, numThetaFreq, numPhiFreq]
                   [0, 1]
                   initVec4 >>= \(plan, vec) ->
                   idft1dGPlan
                     lock
                     plan
                     [numRFreq, numRhoFreq, numThetaFreq, numPhiFreq]
                     [0, 1]
                     vec >>= \(plan, _) ->
                     dft1dGPlan lock plan [numRFreq, numRhoFreq] [0, 1] initVec5 >>= \(plan, vec) ->
                       idft1dGPlan lock plan [numRFreq, numRhoFreq] [0, 1] vec >>= \(plan, _) ->
                         idft1dGPlan
                           lock
                           plan
                           [numThetaFreq, numR2Freq, numR2Freq]
                           [1, 2]
                           initVec6 >>= \(plan, _) ->
                             dft1dGPlan lock plan [numRFreq] [0] initVec7 )
  exportFFTWWisdom (folderPath </> "fftwwisdom.dat")
  return plan

makePlan :: FilePath -> DFTPlan -> Int -> Int -> Int -> Int -> Int -> Int -> Int -> IO DFTPlan
makePlan folderPath initPlan nx ny numR2Freq numThetaFreq numRFreq numPhiFreq numRhoFreq = do
  importFFTWWisdom (folderPath </> "fftwwisdom.dat")
  let lens =
        [ numRFreq * numThetaFreq * numR2Freq ^ 2
        , nx * ny * numThetaFreq * numRFreq
        , numR2Freq ^ 2
        , numThetaFreq * numRFreq * numPhiFreq * numRhoFreq
        , numRhoFreq * numRFreq
        , numThetaFreq * numR2Freq ^ 2
        ]
  (initVec1:initVec2:initVec3:initVec4:initVec5:initVec6:_) <-
   mapConcurrently generateRadomVector lens
  lock <- getFFTWLock
  plan <-
    fst <$>
    (dft1dGPlan
       lock
       initPlan
       [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
       [0, 1, 2, 3]
       initVec1 >>= \(plan, vec) ->
       idft1dGPlan
         lock
         plan
         [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
         [0, 1, 2, 3]
         vec >>= \(plan, _) ->
         dft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] initVec2 >>= \(plan, vec) ->
           idft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] vec -- >>= \(plan, vec) ->
             -- dft1dGPlan
             --   lock
             --   plan
             --   [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
             --   [0]
             --   vec >>= \(plan, vec) ->
             --   idft1dGPlan
             --     lock
             --     plan
             --     [numRFreq, numThetaFreq, numR2Freq, numR2Freq]
             --     [0]
             --     vec
            >>= \(plan, _) ->
             dft1dGPlan lock plan [numR2Freq, numR2Freq] [0, 1] initVec3 >>= \(plan, vec) ->
               idft1dGPlan lock plan [numR2Freq, numR2Freq] [0, 1] vec >>= \(plan, _) ->
                 dft1dGPlan
                   lock
                   plan
                   [numRFreq, numRhoFreq, numThetaFreq, numPhiFreq]
                   [0, 1]
                   initVec4 >>= \(plan, vec) ->
                   idft1dGPlan
                     lock
                     plan
                     [numRFreq, numRhoFreq, numThetaFreq, numPhiFreq]
                     [0, 1]
                     vec >>= \(plan, _) ->
                     dft1dGPlan lock plan [numRFreq, numRhoFreq] [0, 1] initVec5 >>= \(plan, vec) ->
                       idft1dGPlan lock plan [numRFreq, numRhoFreq] [0, 1] vec >>= \(plan, _) ->
                         idft1dGPlan
                           lock
                           plan
                           [numThetaFreq, numR2Freq, numR2Freq]
                           [1, 2]
                           initVec6)
  exportFFTWWisdom (folderPath </> "fftwwisdom.dat")
  return plan

{-# INLINE makePlanFromArray #-}
makePlanFromArray ::
     (R.Source r (Complex Double))
  => DFTPlan
  -> R.Array r DIM4 (Complex Double)
  -> IO DFTPlan
makePlanFromArray initPlan arr = do
  let (Z :. (numRFreq) :. (numThetaFreq) :. (nx) :. (ny)) = extent arr
      initVec = VS.fromList . R.toList $ arr
  lock <- getFFTWLock
  fst <$>
    (dft1dGPlan lock initPlan [nx, ny] [0, 1] initVec >>= \(plan, vec) ->
       idft1dGPlan lock plan [nx, ny] [0, 1] vec >>= \(plan, vec) ->
         dft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] vec >>= \(plan, vec) ->
           idft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] vec)

-- {-# INLINE makePlan #-}
-- makePlan :: FilePath -> DFTPlan -> Int -> Int -> Int -> Int -> Int -> IO DFTPlan
-- makePlan folderPath initPlan nx ny numR2Freq numThetaFreq numRFreq = do
--   importFFTWWisdom (folderPath </> "fftwwisdom.dat")
--   initVec1 <-
--     VS.fromList . L.map (:+ 0) <$> M.replicateM (numR2Freq ^ 2) randomIO
--   initVec2 <-
--     VS.fromList . L.map (:+ 0) <$>
--     M.replicateM (nx * ny * numThetaFreq * numRFreq) randomIO
--   initVec3 <- VS.fromList . L.map (:+ 0) <$> M.replicateM (nx * ny) randomIO
--   lock <- getFFTWLock
--   plan <-
--     fst <$>
--     (dft1dGPlan lock initPlan [numR2Freq, numR2Freq] [0, 1] initVec1 >>= \(plan, vec) ->
--        idft1dGPlan lock plan [numR2Freq, numR2Freq] [0, 1] vec >>= \(plan, _) ->
--          dft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] initVec2 >>= \(plan, vec) ->
--            idft1dGPlan lock plan [numRFreq, numThetaFreq, nx, ny] [0, 1] vec >>= \(plan, _) ->
--              dft1dGPlan lock plan [nx, ny] [0, 1] initVec3 >>= \(plan, vec) ->
--                idft1dGPlan lock plan [nx, ny] [0, 1] vec)
--   exportFFTWWisdom (folderPath </> "fftwwisdom.dat")
--   return plan

ext: hs

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
{-|
Module      : Mandlebrot
Description :
Copyright   : (c) Jason Chen, 2021
License     : MIT
Maintainer  : chezbgone (at) gmail dot com
Stability   : experimental
-}
{-# LANGUAGE NamedFieldPuns #-}
module Mandlebrot where

import Data.Complex
import Data.List (iterate')
import Codec.Picture

import Interpolation

data Resolution = Resolution { width :: Int, height :: Int }

data Bounds = Bounds { top    :: Double
                     , bottom :: Double
                     , left   :: Double
                     , right  :: Double
                     }

data ComputationProps = ComputationProps
    { maxIterations :: Int
    , escapeRadius :: Double
    }

data RenderProps = RenderProps
    { resolution :: Resolution
    , bounds     :: Bounds
    }

data MandlebrotProps = MandlebrotProps
    { computationProps :: ComputationProps
    , renderProps      :: RenderProps
    }

magnitude2 :: RealFloat a => Complex a -> a
magnitude2 (x :+ y) = x*x + y*y

inCardioid :: Complex Double -> Bool
inCardioid (x :+ y) =
    let q = (x - 1/4)^2 + y^2
     in q * (q + x - 1/4) <= 1/4 * y^2

inBulb1 :: Complex Double -> Bool
inBulb1 (x :+ y) = (x + 1)^2 + y^2 <= 1/16

-- https://en.wikipedia.org/wiki/Plotting_algorithms_for_the_Mandelbrot_set#Cardioid_/_bulb_checking
mandlebrotEsc' :: ComputationProps -> Complex Double -> Int
mandlebrotEsc' ComputationProps{ maxIterations, escapeRadius } c =
    let f z = z^2 + c
        orbit = iterate' f 0
     in length $ takeWhile (\p -> magnitude2 p < escapeRadius^2)
              $ take maxIterations orbit

-- optimized
mandlebrotEsc :: ComputationProps -> Complex Double -> Int
mandlebrotEsc cp c = if inCardioid c || inBulb1 c
                         then maxIterations cp
                         else mandlebrotEsc' cp c

xs = (255*) <$> [0, 0.16, 0.42, 0.6425, 0.8575]
rs = [0, 32, 237, 255, 0]
gs = [7, 107, 255, 170, 2]
bs = [100, 203, 255, 0, 0]

rf = interpolate $ piecewiseMonotonicCurve $ zip xs rs
gf = interpolate $ piecewiseMonotonicCurve $ zip xs gs
bf = interpolate $ piecewiseMonotonicCurve $ zip xs bs

intPixel :: Int -> PixelRGB8
intPixel esc =
    let esc' = fromIntegral esc
     in PixelRGB8 (round $ rf esc') (round $ gf esc') (round $ bf esc')

mandlebrotInt :: MandlebrotProps -> DynamicImage
mandlebrotInt mp =
    let cp = computationProps mp
        RenderProps{ resolution=Resolution{ width, height }
                   , bounds=Bounds{ top, bottom, left, right } } = renderProps mp
        dx = (right - left) / fromIntegral width
        dy = (top - bottom) / fromIntegral height
        complex x y = (left + fromIntegral x * dx) :+ (top - fromIntegral y * dy)
     in ImageRGB8 $ generateImage (\x y -> intPixel $ mandlebrotEsc cp $ complex x y) width height

max_stars_repo_name: chezbgone/chaos-renderer

max_issues_repo_name: chezbgone/chaos-renderer

max_forks_repo_name: chezbgone/chaos-renderer


text:
module Main where

import Embadd

import qualified Data.Map.Strict as M (fromList, toList, size, lookup)
import Numeric.LinearAlgebra as LA
import System.Environment

main :: IO ()
main = do
  (vectorFile:args) <- getArgs
  case vectorFile of
    "--help" -> putStrLn "Usage: ./emb.hs FILE"
    "NORM" -> computeNorm (head args)
    _ -> do
      vectorCont <- readFile vectorFile
      let vLine (w:ws) = (w, vector $ map read ws :: Vector Double)
          vectorDict = M.fromList . map (vLine . words) . drop 1 . lines $ vectorCont
      putStrLn "Embedding vizualization tool"
      putStrLn "============================"
      putStrLn ""
      putStrLn "W1 + W2 - W3 -> nearest neighbours of v(W1) + v(W2) - v(W3); you can use any expression with words, '+' and '-'"
      putStrLn "PCA W1 W2 ... Wn  -> pca.png will contain {v(W1), ..., v(Wn)} projected to their principal components"
      readEvalPrintLoop vectorDict

max_stars_repo_name: tomasmcz/embadd

max_issues_repo_name: tomasmcz/embadd

max_forks_repo_name: tomasmcz/embadd


text:
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE Trustworthy #-}
{-# LANGUAGE BangPatterns, DeriveDataTypeable, RecordWildCards #-}
{-# LANGUAGE TypeFamilies #-}

-- |
-- Module      : Gauge.Analysis
-- Copyright   : (c) 2009-2014 Bryan O'Sullivan
--
-- License     : BSD-style
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : GHC
--
-- Analysis code for benchmarks.

-- XXX Do we really want to expose this module to users? This is all internal.
-- Most of the exports are not even used by Gauge itself outside of this
-- module.

module Gauge.Analysis
    ( Outliers(..)
    , OutlierEffect(..)
    , OutlierVariance(..)
    , Report
    , SampleAnalysis(..)
    , analyseSample
    , scale
    , analyseBenchmark
    , analyseMean
    , countOutliers
    , classifyOutliers
    , noteOutliers
    , outlierVariance
    , regress
    , benchmark'
    , benchmarkWith'
    ) where

import Control.Applicative
import Control.Arrow (second)
import Control.DeepSeq (NFData(rnf))
import Control.Monad (forM_, when)
import Gauge.Benchmark (Benchmark (..), BenchmarkAnalysis(..), Benchmarkable, runBenchmark)
import Gauge.IO.Printf (note, printError, prolix, rewindClearLine)
import Gauge.Main.Options (defaultConfig, Config(..), Verbosity (..),
                           DisplayMode (..))
import Gauge.Measurement (Measured(measTime), secs, rescale, measureKeys,
                          measureAccessors_, validateAccessors, renderNames)
import Gauge.Monad (Gauge, askConfig, gaugeIO, Crit(..), askCrit, withConfig)
import Gauge.Format
import qualified Gauge.CSV as CSV
import Data.Data (Data, Typeable)
import Data.Int (Int64)
import Data.IORef (IORef, readIORef, writeIORef)
import Gauge.ListMap (Map)
import Data.Maybe (fromJust, isJust)
import Data.Traversable
import GHC.Generics (Generic)
import Statistics.Function (sort)
import Statistics.Quantile (weightedAvg, Sorted(..))
import Statistics.Regression (bootstrapRegress, olsRegress)
import Statistics.Resampling (Estimator(..),resample)
import Statistics.Sample (mean)
import Statistics.Sample.KernelDensity (kde)
import Statistics.Types (Sample, Estimate(..),ConfInt(..),confidenceInterval
                        ,cl95,confidenceLevel)
import System.Random.MWC (GenIO, createSystemRandom)
import Text.Printf (printf)
import qualified Gauge.ListMap as Map
import qualified Data.Vector as V
import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import qualified Statistics.Resampling.Bootstrap as B
import qualified Statistics.Types                as B
import qualified Statistics.Types as St
import Prelude hiding (sequence, mapM)

-- | Outliers from sample data, calculated using the boxplot
-- technique.
data Outliers = Outliers {
      samplesSeen :: !Int64
    , lowSevere   :: !Int64
    -- ^ More than 3 times the interquartile range (IQR) below the
    -- first quartile.
    , lowMild     :: !Int64
    -- ^ Between 1.5 and 3 times the IQR below the first quartile.
    , highMild    :: !Int64
    -- ^ Between 1.5 and 3 times the IQR above the third quartile.
    , highSevere  :: !Int64
    -- ^ More than 3 times the IQR above the third quartile.
    } deriving (Eq, Show, Typeable, Data, Generic)

instance NFData Outliers

-- | A description of the extent to which outliers in the sample data
-- affect the sample mean and standard deviation.
data OutlierEffect = Unaffected -- ^ Less than 1% effect.
                   | Slight     -- ^ Between 1% and 10%.
                   | Moderate   -- ^ Between 10% and 50%.
                   | Severe     -- ^ Above 50% (i.e. measurements
                                -- are useless).
                     deriving (Eq, Ord, Show, Typeable, Data, Generic)

instance NFData OutlierEffect

outliersEmpty :: Outliers
outliersEmpty = Outliers 0 0 0 0 0

addOutliers :: Outliers -> Outliers -> Outliers
addOutliers (Outliers s a b c d) (Outliers t w x y z) =
    Outliers (s+t) (a+w) (b+x) (c+y) (d+z)
{-# INLINE addOutliers #-}

-- | Analysis of the extent to which outliers in a sample affect its
-- standard deviation (and to some extent, its mean).
data OutlierVariance = OutlierVariance {
      ovEffect   :: OutlierEffect
    -- ^ Qualitative description of effect.
    , ovDesc     :: String
    -- ^ Brief textual description of effect.
    , ovFraction :: Double
    -- ^ Quantitative description of effect (a fraction between 0 and 1).
    } deriving (Eq, Show, Typeable, Data, Generic)

instance NFData OutlierVariance where
    rnf OutlierVariance{..} = rnf ovEffect `seq` rnf ovDesc `seq` rnf ovFraction

-- | Results of a linear regression.
data Regression = Regression {
    regResponder  :: String
    -- ^ Name of the responding variable.
  , regCoeffs     :: Map String (St.Estimate St.ConfInt Double)
    -- ^ Map from name to value of predictor coefficients.
  , regRSquare    :: St.Estimate St.ConfInt Double
    -- ^ R² goodness-of-fit estimate.
  } deriving (Eq, Show, Typeable, Generic)

instance NFData Regression where
    rnf Regression{..} =
      rnf regResponder `seq` rnf regCoeffs `seq` rnf regRSquare

-- | Result of a bootstrap analysis of a non-parametric sample.
data SampleAnalysis = SampleAnalysis {
      anRegress    :: [Regression]
      -- ^ Estimates calculated via linear regression.
    , anMean       :: St.Estimate St.ConfInt Double
      -- ^ Estimated mean.
    , anStdDev     :: St.Estimate St.ConfInt Double
      -- ^ Estimated standard deviation.
    , anOutlierVar :: OutlierVariance
      -- ^ Description of the effects of outliers on the estimated
      -- variance.
    } deriving (Eq, Show, Typeable, Generic)

instance NFData SampleAnalysis where
    rnf SampleAnalysis{..} =
        rnf anRegress `seq` rnf anMean `seq`
        rnf anStdDev `seq` rnf anOutlierVar

-- | Data for a KDE chart of performance.
data KDE = KDE {
      kdeType   :: String
    , kdeValues :: U.Vector Double
    , kdePDF    :: U.Vector Double
    } deriving (Eq, Show, Typeable, Data, Generic)

instance NFData KDE where
    rnf KDE{..} = rnf kdeType `seq` rnf kdeValues `seq` rnf kdePDF

-- | Report of a sample analysis.
data Report = Report {
      reportName     :: String
      -- ^ The name of this report.
    , reportKeys     :: [String]
      -- ^ See 'measureKeys'.
    , reportMeasured :: V.Vector Measured
      -- ^ Raw measurements.
    , reportAnalysis :: SampleAnalysis
      -- ^ Report analysis.
    , reportOutliers :: Outliers
      -- ^ Analysis of outliers.
    , reportKDEs     :: [KDE]
      -- ^ Data for a KDE of times.
    } deriving (Eq, Show, Typeable, Generic)

instance NFData Report where
    rnf Report{..} =
      rnf reportName `seq` rnf reportKeys `seq`
      rnf reportMeasured `seq` rnf reportAnalysis `seq` rnf reportOutliers `seq`
      rnf reportKDEs

-- | Classify outliers in a data set, using the boxplot technique.
classifyOutliers :: Sample -> Outliers
classifyOutliers sa = U.foldl' ((. outlier) . addOutliers) outliersEmpty ssa
    where outlier e = Outliers {
                        samplesSeen = 1
                      , lowSevere = if e <= loS && e < hiM then 1 else 0
                      , lowMild = if e > loS && e <= loM then 1 else 0
                      , highMild = if e >= hiM && e < hiS then 1 else 0
                      , highSevere = if e >= hiS && e > loM then 1 else 0
                      }
          !loS = q1 - (iqr * 3)
          !loM = q1 - (iqr * 1.5)
          !hiM = q3 + (iqr * 1.5)
          !hiS = q3 + (iqr * 3)
          q1   = weightedAvg 1 4 (Sorted ssa)
          q3   = weightedAvg 3 4 (Sorted ssa)
          ssa  = sort sa
          iqr  = q3 - q1

-- | Compute the extent to which outliers in the sample data affect
-- the sample mean and standard deviation.
outlierVariance
  :: B.Estimate B.ConfInt Double -- ^ Bootstrap estimate of sample mean.
  -> B.Estimate B.ConfInt Double -- ^ Bootstrap estimate of sample
                                 --   standard deviation.
  -> Double                      -- ^ Number of original iterations.
  -> OutlierVariance
outlierVariance µ σ a = OutlierVariance effect desc varOutMin
  where
    ( effect, desc ) | varOutMin < 0.01 = (Unaffected, "no")
                     | varOutMin < 0.1  = (Slight,     "slight")
                     | varOutMin < 0.5  = (Moderate,   "moderate")
                     | otherwise        = (Severe,     "severe")
    varOutMin = (minBy varOut 1 (minBy cMax 0 µgMin)) / σb2
    varOut c  = (ac / a) * (σb2 - ac * σg2) where ac = a - c
    σb        = B.estPoint σ
    µa        = B.estPoint µ / a
    µgMin     = µa / 2
    σg        = min (µgMin / 4) (σb / sqrt a)
    σg2       = σg * σg
    σb2       = σb * σb
    minBy f q r = min (f q) (f r)
    cMax x    = fromIntegral (floor (-2 * k0 / (k1 + sqrt det)) :: Int)
      where
        k1    = σb2 - a * σg2 + ad
        k0    = -a * ad
        ad    = a * d
        d     = k * k where k = µa - x
        det   = k1 * k1 - 4 * σg2 * k0

-- | Count the total number of outliers in a sample.
countOutliers :: Outliers -> Int64
countOutliers (Outliers _ a b c d) = a + b + c + d
{-# INLINE countOutliers #-}

-- | Display the mean of a 'Sample', and characterise the outliers
-- present in the sample.
analyseMean :: Sample
            -> Int              -- ^ Number of iterations used to
                                -- compute the sample.
            -> Gauge Double
analyseMean a iters = do
  let µ = mean a
  _ <- note "mean is %s (%d iterations)\n" (secs µ) iters
  noteOutliers . classifyOutliers $ a
  return µ

-- | Multiply the 'Estimate's in an analysis by the given value, using
-- 'B.scale'.
scale :: Double                 -- ^ Value to multiply by.
      -> SampleAnalysis -> SampleAnalysis
scale f s@SampleAnalysis{..} = s {
                                 anMean = B.scale f anMean
                               , anStdDev = B.scale f anStdDev
                               }

-- | Return a random number generator, creating one if necessary.
--
-- This is not currently thread-safe, but in a harmless way (we might
-- call 'createSystemRandom' more than once if multiple threads race).
getGen :: Gauge GenIO
getGen = memoise gen createSystemRandom

getMeasurement :: (U.Unbox a) => (Measured -> a) -> V.Vector Measured -> U.Vector a
getMeasurement f v = U.convert . V.map f $ v

-- | Memoise the result of an 'IO' action.
--
-- This is not currently thread-safe, but hopefully in a harmless way.
-- We might call the given action more than once if multiple threads
-- race, so our caller's job is to write actions that can be run
-- multiple times safely.
memoise :: (Crit -> IORef (Maybe a)) -> IO a -> Gauge a
memoise ref generate = do
  r <- ref <$> askCrit
  gaugeIO $ do
    mv <- readIORef r
    case mv of
      Just rv -> return rv
      Nothing -> do
        rv <- generate
        writeIORef r (Just rv)
        return rv

toCL :: Maybe Double -> B.CL Double
toCL a =
    case a of
        Nothing -> B.cl95
        Just x -> B.mkCL x

-- | Perform an analysis of a measurement.
analyseSample :: String            -- ^ Experiment name.
              -> V.Vector Measured -- ^ Sample data.
              -> Gauge (Either String Report)
analyseSample name meas = do
  Config{..} <- askConfig
  let ests      = [Mean,StdDev]
      stime     = getMeasurement (measTime . rescale) $ meas
      n         = G.length meas
  _ <- prolix "bootstrapping with %d samples\n" n

  gen <- getGen
  ers <- (sequence <$>) . mapM (\(ps,r) -> regress gen ps r meas) $ ((["iters"],"time"):regressions)
  case ers of
    Left err -> pure $ Left err
    Right rs -> do
      resamps <- gaugeIO $ resample gen ests resamples stime
      let [estMean,estStdDev] = B.bootstrapBCA (toCL confInterval) stime
                                               resamps
          ov = outlierVariance estMean estStdDev (fromIntegral n)
          an = SampleAnalysis
                 { anRegress    = rs
                 , anMean       = estMean
                 , anStdDev     = estStdDev
                 , anOutlierVar = ov
                 }
      return $ Right $ Report
        { reportName     = name
        , reportKeys     = measureKeys
        , reportMeasured = meas
        , reportAnalysis = an
        , reportOutliers = classifyOutliers stime
        , reportKDEs     = [uncurry (KDE "time") (kde 128 stime)]
        }

-- | Regress the given predictors against the responder.
--
-- Errors may be returned under various circumstances, such as invalid
-- names or lack of needed data.
--
-- See 'olsRegress' for details of the regression performed.
regress :: GenIO
        -> [String]             -- ^ Predictor names.
        -> String               -- ^ Responder name.
        -> V.Vector Measured
        -> Gauge (Either String Regression)
regress gen predNames respName meas
    | G.null meas = pure $ Left "no measurements"
    | otherwise   = case validateAccessors predNames respName of
        Left err   -> pure $ Left err
        Right accs -> do
            let unmeasured = [n | (n, Nothing) <- map (second ($ G.head meas)) accs]
            if not (null unmeasured)
                then pure $ Left $ "no data available for " ++ renderNames unmeasured
                else do
                    let (r:ps) = map ((`getMeasurement` meas) . (fromJust .) . snd) accs
                    Config{..} <- askConfig
                    (coeffs,r2) <- gaugeIO $
                        bootstrapRegress gen resamples (toCL confInterval)
                                         olsRegress ps r
                    pure $ Right $ Regression
                        { regResponder = respName
                        , regCoeffs    = Map.fromList (zip (predNames ++ ["y"]) (G.toList coeffs))
                        , regRSquare   = r2
                        }

-- | Display a report of the 'Outliers' present in a 'Sample'.
noteOutliers :: Outliers -> Gauge ()
noteOutliers o = do
  let frac n = (100::Double) * fromIntegral n / fromIntegral (samplesSeen o)
      check :: Int64 -> Double -> String -> Gauge ()
      check k t d = when (frac k > t) $
                    note "  %d (%.1g%%) %s\n" k (frac k) d
      outCount = countOutliers o
  when (outCount > 0) $ do
    _ <- note "found %d outliers among %d samples (%.1g%%)\n"
         outCount (samplesSeen o) (frac outCount)
    check (lowSevere o) 0 "low severe"
    check (lowMild o) 1 "low mild"
    check (highMild o) 1 "high mild"
    check (highSevere o) 0 "high severe"

-- | Analyse a single benchmark.
analyseBenchmark :: String -> V.Vector Measured -> Gauge Report
analyseBenchmark desc meas = do
  Config{..} <- askConfig
  _ <- prolix "%sanalysing with %d resamples\n" rewindClearLine resamples
  -- XXX handle meas being empty
  erp <- analyseSample desc meas
  case erp of
    Left err -> printError "*** Error: %s\n" err
    Right rpt@Report{..} -> do
        let SampleAnalysis{..} = reportAnalysis
            OutlierVariance{..} = anOutlierVar
            wibble = printOverallEffect ovEffect
            (builtin, others) = splitAt 1 anRegress

        gaugeIO $ CSV.write csvFile $ CSV.Row
            [ CSV.string desc
            , CSV.float (estPoint anMean)
            , CSV.float (fst $ confidenceInterval anMean)
            , CSV.float (snd $ confidenceInterval anMean)
            , CSV.float (estPoint anStdDev)
            , CSV.float (fst $ confidenceInterval anStdDev)
            , CSV.float (snd $ confidenceInterval anStdDev)
            ]

        case displayMode of
            StatsTable -> do
              _ <- note "%sbenchmarked %s%s%s\n" rewindClearLine green desc reset
              let r2 n = printf "%.3f R\178" n
              forM_ builtin $ \Regression{..} ->
                case Map.lookup "iters" regCoeffs of
                  Nothing -> return ()
                  Just t  -> bs secs "time" t >> bs r2 "" regRSquare
              bs secs "mean" anMean
              bs secs "std dev" anStdDev
              forM_ others $ \Regression{..} -> do
                _ <- bs r2 (regResponder ++ ":") regRSquare
                forM_ (Map.toList regCoeffs) $ \(prd,val) ->
                  bs (printf "%.3g") ("  " ++ prd) val
              when (verbosity == Verbose || (ovEffect > Slight && verbosity > Quiet)) $ do
                when (verbosity == Verbose) $ noteOutliers reportOutliers
                _ <- note "variance introduced by outliers: %d%% (%s)\n"
                     (round (ovFraction * 100) :: Int) wibble
                return ()

              _ <- traverse
                    (\(k, (a, s, _)) -> reportStat Verbose a s k)
                    measureAccessors_
              _ <- note "\n"
              pure ()
            Condensed -> do
              _ <- note "%s%s%-40s%s " rewindClearLine green desc reset
              bsSmall secs "mean" anMean
              bsSmall secs "( +-" anStdDev
              _ <- note ")\n"
              pure ()

        return rpt
      where bs :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge ()
            bs f metric e@Estimate{..} =
              note "%s%-20s%s %s%-10s%s (%s .. %s%s)\n" yellow metric reset
                   red (f estPoint) reset
                   (f $ fst $ confidenceInterval e) (f $ snd $ confidenceInterval e)
                   (let cl = confIntCL estError
                        str | cl == cl95 = ""
                            | otherwise  = printf ", ci %.3f" (confidenceLevel cl)
                    in str
                   )
            bsSmall :: (Double -> String) -> String -> Estimate ConfInt Double -> Gauge ()
            bsSmall f metric Estimate{..} =
              note "%s %s%-10s%s" metric red (f estPoint) reset

            reportStat :: Verbosity
                       -> (Measured -> Maybe Double)
                       -> (Double -> String)
                       -> String -> Gauge ()
            reportStat lvl accessor sh msg = do
              let v0 = V.map (accessor . rescale) meas
              -- Print average only if all data points are present
              when (V.all isJust v0) $ do
                  let v = V.map fromJust v0
                      total = V.sum v
                      len = V.length v
                      avg = total / (fromIntegral len)
                  when (verbosity >= lvl && avg > 0.0) $ do
                    note "%-20s %-10s (%s .. %s)\n" msg (sh avg)
                      (sh (V.minimum v)) (sh (V.maximum v))

printOverallEffect :: OutlierEffect -> String
printOverallEffect Unaffected = "unaffected"
printOverallEffect Slight     = "slightly inflated"
printOverallEffect Moderate   = "moderately inflated"
printOverallEffect Severe     = "severely inflated"

-- XXX The original type of these types returned 'Report' type. But the
-- implementation was wrong as it was not running any environment settings on
-- the way to the benchmark. Now, we have used the correct function to do that
-- but unfortunately that function returns void. That can be fixed though if it
-- is important.

-- | Run a benchmark interactively and analyse its performance.
benchmarkWith' :: Config -> Benchmarkable -> IO ()
benchmarkWith' cfg bm = withConfig cfg $
    runBenchmark (const True) (Benchmark "function" bm) (BenchmarkNormal analyseBenchmark)

-- | Run a benchmark interactively and analyse its performanc.
benchmark' :: Benchmarkable -> IO ()
benchmark' = benchmarkWith' defaultConfig

max_stars_repo_name: runeksvendsen/hs-gauge

max_issues_repo_name: runeksvendsen/hs-gauge

max_forks_repo_name: runeksvendsen/hs-gauge


text:
{-# OPTIONS_GHC -Wno-incomplete-patterns #-}

module NVectorSpec where

import Data.Complex

import Test.QuickCheck
import Test.QuickCheck.Poly

import Category
import Comonad
import Functor
import NVector
import Unboxed

type N = 10

prop_NVector_Functor_id :: FnProp (NVector N A)
prop_NVector_Functor_id = checkFnEqual law_Functor_id

prop_NVector_Functor_comp :: Fun B C -> Fun A B -> FnProp (NVector N A)
prop_NVector_Functor_comp f g = checkFnEqual (law_Functor_comp f g)

prop_NVector_Semicomonad_comm ::
    Fun (NVector N B) C -> Fun (NVector N A) B -> FnProp (NVector N A)
prop_NVector_Semicomonad_comm f g = checkFnEqual (law_Semicomonad_comm f g)

prop_NVector_Semicomonad1_comm ::
    Fun (NVector N B) C -> Fun (NVector N A) B -> FnProp (NVector N A)
prop_NVector_Semicomonad1_comm f g = checkFnEqual (law_Semicomonad1_comm f g)

prop_NVector_Semicomonad1_comm' ::
    Fun (NVector N B) C -> Fun (NVector N A) B -> FnProp (NVector N A)
prop_NVector_Semicomonad1_comm' (Fn f) (Fn g) =
    checkFnEqual (law_Semicomonad1_comm' f g)

type UA = Int
type UB = Double
type UC = Complex Double

prop_NUVector_Functor_id :: FnProp (NUVector N UA)
prop_NUVector_Functor_id = checkFnEqual law_Functor_id

prop_NUVector_Functor_comp ::
    (UB -#> UC) -> (UA -#> UB) -> FnProp (NUVector N UA)
prop_NUVector_Functor_comp f g = checkFnEqual (law_Functor_comp f g)

prop_NUVector_Semicomonad1_comm' ::
    Fun (NUVector N UB) UC -> Fun (NUVector N UA) UB ->
    FnProp (NUVector N UA)
prop_NUVector_Semicomonad1_comm' (Fn f) (Fn g) =
    checkFnEqual (law_Semicomonad1_comm' f g)


text:
{-# language ScopedTypeVariables #-}
{-# language TypeFamilies #-}
module RealMain where

import Codec.Picture
import Control.Monad.State
import Data.Word
import Data.Complex
import System.Exit (exitFailure, exitSuccess)
import System.Environment (getArgs)
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Generic as G

import Paths_proyecto_lena
import ArrayPoking
import Algoritmo

realMain :: IO ()
realMain =
  do rate : divisions : reconstruct : _ <- getArgs
     img@(Image width height _) <- loadImage lena_path
     let red   = V.map fromIntegral (extractColor 0 img)
         green = V.map fromIntegral (extractColor 1 img)
         blue  = V.map fromIntegral (extractColor 2 img)

         calidad    = read rate
         divisiones = read divisions
         recons     = read reconstruct
         factor     = divisiones - recons

         redNew   = waveletAlgo calidad divisiones recons red
         greenNew = waveletAlgo calidad divisiones recons green
         blueNew  = waveletAlgo calidad divisiones recons blue

     writePng "./lena_alt.png"
              (nuevaImg factor width height redNew greenNew blueNew)

realMain2 :: IO ()
realMain2 =
  do rate : divisions : reconstruct : _ <- getArgs
     img@(Image width height _) <- loadImage lena_path
     let red   = V.map fromIntegral (extractColor 0 img)

         calidad    = read rate
         divisiones = read divisions
         recons     = read reconstruct
         factor     = divisiones - recons

         redNew   = waveletAlgo calidad divisiones recons red

     writePng "./lena_alt.png"
       -- Solo afecta el largo pues es 1D
       (Image (width `div` (2^(factor - 1))) (height)
              (V.map (round . realPart) redNew) :: Image Pixel8)

loadImage :: FilePath -> IO (Image PixelRGB8)
loadImage path =
  do lena <- readImage path
     img  <- either (\s -> putStrLn s >> exitFailure) return lena
     case img of
       ImageRGB8 img_lena -> return img_lena
       _                  -> exitFailure

loadImageBebe :: FilePath -> IO (Image Pixel8)
loadImageBebe path =
  do lena <- readImage path
     img  <- either (\s -> putStrLn s >> exitFailure) return lena
     case img of
       ImageY8   img_bebe -> return img_bebe
       _                  -> exitFailure

nuevaImg :: (Complex Double ~ a)
         => Int -> Int -> Int
         -> V.Vector a -> V.Vector a -> V.Vector a
         -> Image PixelRGB8
nuevaImg factor width height red green blue =
  let combined = combineColors red green blue
   in Image (width `div` 2^(factor - 1)) height
      . G.map (round . realPart) $ combined

lena_path, tiny_path, bebe_path :: FilePath
lena_path = "/home/slack/UTFSM/2016-2/trabajo-roldan-peypo/proyecto-lena/images/Lenna.png"
tiny_path = "/home/slack/UTFSM/2016-2/trabajo-roldan-peypo/proyecto-lena/images/tiny.png"
bebe_path = "/home/slack/UTFSM/2016-2/trabajo-roldan-peypo/proyecto-lena/images/bebe.png"

max_stars_repo_name: RubenAstudillo/proyecto-lena

max_issues_repo_name: RubenAstudillo/proyecto-lena

max_forks_repo_name: RubenAstudillo/proyecto-lena


text:
module PlotLocalEigenVector where

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

main = do
  args@(numPointStr:numOrientationStr:numScaleStr:thetaSigmaStr:scaleSigmaStr:maxScaleStr:taoStr:numTrailStr:maxTrailStr:thetaFreqsStr:scaleFreqsStr:histFileName:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      numScale = read numScaleStr :: Int
      thetaSigma = read thetaSigmaStr :: Double
      scaleSigma = read scaleSigmaStr :: Double
      maxScale = read maxScaleStr :: Double
      tao = read taoStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      scaleFreq = read scaleFreqsStr :: Double
      scaleFreqs = [-scaleFreq .. scaleFreq]
      numThread = read numThreadStr :: Int
      folderPath = "output/test/PlotLocalEigenVector"
      histFilePath = folderPath </> histFileName
      imageFolder =
        printf
          "%d_%.2f_%.2f_%.0f_%.0f_%.0f_%.0f"
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          thetaFreq
          scaleFreq
  removePathForcibly (folderPath </> imageFolder)
  createDirectoryIfMissing True (folderPath </> imageFolder)
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z2T0S0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          thetaSigma
          scaleSigma
          maxScale
          tao
          thetaFreqs
          thetaFreqs
          scaleFreqs
          scaleFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length scaleFreqs
             , L.length thetaFreqs
             , L.length scaleFreqs
             , L.length thetaFreqs
             ]
             0)
  let localEigenVector =
        computeLocalEigenVector
          (ParallelParams numThread 0)
          (pinwheelHollow 8)
          (cutoff 64 radialArr)
          numPoint
          numPoint
          maxScale
          thetaFreqs
          scaleFreqs
      rotationDeg = 90
      rotatedLocalEigenVector =
        R.traverse2
          localEigenVector
          (fromListUnboxed (Z :. (L.length thetaFreqs)) thetaFreqs)
          const $ \f fFreq idx@(Z :. k :. l :. i :. j) ->
          f idx *
          ((exp (0 :+ (-1) * (fFreq (Z :. k)) * rotationDeg / 180 * pi)) --  +
           -- (exp (0 :+ (-1) * (fFreq (Z :. k)) * (rotationDeg + 90) / 180 * pi)) +
           -- (exp (0 :+ (-1) * (fFreq (Z :. k)) * (rotationDeg + 180) / 180 * pi)) +
           -- (exp (0 :+ (-1) * (fFreq (Z :. k)) * (rotationDeg + 270) / 180 * pi))
           )
  MP.mapM_
    (\(tf, sf) ->
       plotImageRepaComplex
         (folderPath </> imageFolder </> (printf "%d_%d.png" tf sf)) .
       ImageRepa 8 .
       computeS .
       R.extend (Z :. (1 :: Int) :. All :. All) . R.slice localEigenVector $
       (Z :. tf :. sf :. All :. All))
    [ (tf, sf)
    | tf <- [0 .. L.length thetaFreqs - 1]
    , sf <- [0 .. L.length scaleFreqs - 1]
    ]
  plotImageRepaComplex (folderPath </> imageFolder </> "sumFreq.png") .
    ImageRepa 8 .
    computeS .
    R.extend (Z :. (1 :: Int) :. All :. All) .
    R.sumS . R.sumS . rotate4D . rotate4D $
    localEigenVector
  plotImageRepaComplex (folderPath </> imageFolder </> "sumFreqRotation.png") .
    ImageRepa 8 .
    computeS .
    R.extend (Z :. (1 :: Int) :. All :. All) .
    R.sumS . R.sumS . rotate4D . rotate4D $
    rotatedLocalEigenVector
  let localEigenVectorR2S1RP =
        r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
        localEigenVector
      rotatedLocalEigenVectorR2S1RP =
        r2z2Tor2s1rp numOrientation thetaFreqs numScale scaleFreqs $
        rotatedLocalEigenVector
  plotImageRepaComplex (folderPath </> imageFolder </> "sum.png") .
    ImageRepa 8 .
    computeS .
    R.extend (Z :. (1 :: Int) :. All :. All) .
    R.sumS . R.sumS . rotate4D . rotate4D $
    localEigenVectorR2S1RP
  plotImageRepaComplex (folderPath </> imageFolder </> "sumRotation.png") .
    ImageRepa 8 .
    computeS .
    R.extend (Z :. (1 :: Int) :. All :. All) .
    R.sumS . R.sumS . rotate4D . rotate4D $
    localEigenVectorR2S1RP
  MP.mapM_
    (\i ->
       plotImageRepaComplex
         (folderPath </> imageFolder </> printf "Theta%03d.png" i) .
       ImageRepa 8 .
       computeS .
       R.extend (Z :. (1 :: Int) :. All :. All) . R.slice localEigenVectorR2S1RP $
       (Z :. i :. (0 :: Int) :. All :. All))
    [0 .. numOrientation - 1]


text:
{-# LANGUAGE DataKinds  #-}
{-# LANGUAGE RankNTypes #-}

module Main where

import qualified Data.Vector.Sized as V
import           GHC.TypeLits (KnownNat)
import           Graphics.Gloss
import           Numeric.Hamilton
import           Numeric.LinearAlgebra.Static
import           Pendulum.Model
import           Pendulum.View

toVector :: forall a m. KnownNat m => [a] -> V.Vector m a
toVector = (\(Just x) -> x) . V.fromList

pdConfig :: Config 4
pdConfig =
  let th = 0.5 * pi
   in Cfg (vector (replicate 4 th)) (vector (replicate 4 0))

pendulum :: Pendulum 4
pendulum = Pendulum
  { _lengthes = toVector (replicate 4 0.8)
  , _masses   = toVector (replicate 4 0.8)
  , _config   = pdConfig
  }

pdSystem :: System 8 4
pdSystem = mkPendlumSystem pendulum

main :: IO ()
main = do
  let display = InWindow "Multi Pendulum Simulation" (800, 800) (10, 10)
  simulate display white 30 pendulum drawPendulum (\_ dt pd -> stepPendulum (realToFrac dt) pdSystem pd)


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

module Polynomials where

import           Numeric.LinearAlgebra
import           System.Random hiding (rand)
import           Control.Monad.IO.Class
import           Control.Monad
import           Control.Exception
import           Data.Maybe
import qualified Data.Vector as V
import qualified Data.Array.Accelerate as A
import           Genetic
import           Formulas

import Debug.Trace

data Article = Article
            { _articlePowers  :: Powers
            , _articleInputs  :: Inplist
            , _articleOutputs :: Y
            } deriving (Eq, Ord)

instance Show Article where
      show = show . _articlePowers

instance Individual Article where
   -- fitness :: a -> IO (Maybe Double)
      fitness (Article ps is os) = do
            let x' = calcX is ps
                rows = length $ toLists x'
                cols = length . head $ toLists x'
                y = os
            noise_ <- rand rows cols
            let   noise = noise_ * 100.0
                  x = if det (tr x' <> x') == 0
                      then x' + noise
                      else x'
            -- bug fix for BLAS matrix rank on big numbers
            return $ if det (tr x <> x) == 0
                     then Nothing
                     else Just $ cost (normEq x y) x y

   -- mutateIndividual :: a -> IO a
      mutateIndividual (Article pows _a _b) = do
            r <- randomIO
            if mutationProbability < r
            then inner pows
            else return $ Article pows _a _b
            where
                  mutationProbability :: Double
                  mutationProbability = 0.1

                  inner :: Powers -> IO Article
                  inner x = do
                        blockIndex <- randomRIO (0, V.length x - 1)
                        powerIndex <- randomRIO (0, V.length (x V.! 0) - 1)
                        newPower   <- randomRIO (-4, 4)

                        let updated = updateVV x ( blockIndex, powerIndex) newPower
                        return $ Article updated _a _b

   -- crossover :: (a, a) -> IO a
      crossover ( Article pows1 _a _b
                , Article pows2 _ _
                ) = do
            index <- randomRIO (0, V.length pows1)
            let leftPart  = V.take index pows1
            let rightPart = V.drop index pows2
            return $ Article (V.concat [leftPart, rightPart]) _a _b

updateVV :: V.Vector (V.Vector a) -- what to update
         -> (Int, Int)            -- where
         -> a                     -- what to put there
         -> V.Vector (V.Vector a) -- result
updateVV vvs (xIdx, yIdx) value =
      vvs V.// [(xIdx, updatedV)]
      where
            updatedV  = originalV V.// [(yIdx, value)]
            originalV = vvs V.! xIdx


text:
{-# OPTIONS_GHC -XTemplateHaskell #-}

module Euler.Problem011Test (suite) where

import Numeric.LinearAlgebra.Data (Matrix, loadMatrix)
import Test.Tasty (TestTree)
import Test.Tasty.HUnit
import Test.Tasty.TH (testGroupGenerator)

import Euler.Problem011

suite :: TestTree
suite = $(testGroupGenerator)

matrix :: IO (Matrix Double)
matrix = loadMatrix "input/011.txt"

case_count_of_vectors_in_input :: Assertion
case_count_of_vectors_in_input =
  matrix >>= \m -> 10 @=? (length . vectors . head $ submatrices 4 m)

case_count_of_4x4_submatrices_in_input :: Assertion
case_count_of_4x4_submatrices_in_input =
  matrix >>= \m -> 17 * 17 @=? (length $ submatrices 4 m)

max_stars_repo_name: whittle/euler

max_issues_repo_name: whittle/euler

max_forks_repo_name: whittle/euler


text:
{-# LANGUAGE BangPatterns        #-}
{-# LANGUAGE CPP                 #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TupleSections       #-}
{-# LANGUAGE TypeFamilies        #-}
{-# LANGUAGE TypeOperators       #-}

import           Control.Applicative
import           Control.DeepSeq
import           Control.Monad
import           Control.Monad.Random
import           Control.Monad.Trans.Except

import qualified Data.Attoparsec.Text         as A
import           Data.List                    (foldl')
import qualified Data.Text                    as T
import qualified Data.Text.IO                 as T
import qualified Data.Vector.Storable         as V

import           Numeric.LinearAlgebra        (maxIndex)
import qualified Numeric.LinearAlgebra.Static as SA

import           Options.Applicative

import           Grenade
import           Grenade.Utils.OneHot

--
-- Note: Input files can be downloaded at https://www.kaggle.com/scolianni/mnistasjpg
--

-- It's logistic regression!
--
-- This network is used to show how we can embed a Network as a layer in the larger MNIST
-- type.
type FL i o =
  Network
    '[ FullyConnected i o, Logit ]
    '[ 'D1 i, 'D1 o, 'D1 o ]

-- The definition of our convolutional neural network.
-- In the type signature, we have a type level list of shapes which are passed between the layers.
-- One can see that the images we are inputing are two dimensional with 28 * 28 pixels.

-- It's important to keep the type signatures, as there's many layers which can "squeeze" into the gaps
-- between the shapes, so inference can't do it all for us.

-- With the mnist data from Kaggle normalised to doubles between 0 and 1, learning rate of 0.01 and 15 iterations,
-- this network should get down to about a 1.3% error rate.
--
-- /NOTE:/ This model is actually too complex for MNIST, and one should use the type given in the readme instead.
--         This one is just here to demonstrate Inception layers in use.
--
type MNIST =
  Network
    '[ Reshape,
       Concat ('D3 28 28 1) Trivial ('D3 28 28 14) (InceptionMini 28 28 1 5 9),
       Pooling 2 2 2 2, Relu,
       Concat ('D3 14 14 3) (Convolution 15 3 1 1 1 1) ('D3 14 14 15) (InceptionMini 14 14 15 5 10), Crop 1 1 1 1, Pooling 3 3 3 3, Relu,
       Reshape, FL 288 80, FL 80 10 ]
    '[ 'D2 28 28, 'D3 28 28 1,
       'D3 28 28 15, 'D3 14 14 15, 'D3 14 14 15, 'D3 14 14 18,
       'D3 12 12 18, 'D3 4 4 18, 'D3 4 4 18,
       'D1 288, 'D1 80, 'D1 10 ]

randomMnist :: IO MNIST
randomMnist = randomNetwork

convTest :: Int -> FilePath -> FilePath -> Optimizer opt -> ExceptT String IO ()
convTest iterations trainFile validateFile opt = do
  net0         <- lift randomMnist
  trainData    <- readMNIST trainFile
  validateData <- readMNIST validateFile
  lift $ foldM_ (runIteration trainData validateData) net0 [1..iterations]

    where
  trainEach !opt' !network (!i, !o) = force (train opt' network i o)

  runIteration !trainRows !validateRows !net !i = do
    putStrLn $ "Number of training rows: " ++ show (length trainRows)
    let !trained' = foldl' (trainEach (sgdUpdateLearningParamters opt)) net trainRows
    let !res      = fmap (\(rowP,rowL) -> (rowL,) $ runNet trained' rowP) validateRows
    let !res'     = fmap (\(S1D label, S1D prediction) -> (maxIndex (SA.extract label), maxIndex (SA.extract prediction))) res
    print trained'
    putStrLn $ "Iteration " ++ show i ++ ": " ++ show (length (filter ((==) <$> fst <*> snd) res')) ++ " of " ++ show (length res')
    return trained'
  sgdUpdateLearningParamters :: Optimizer opt -> Optimizer opt
  sgdUpdateLearningParamters (OptSGD rate mom reg) = OptSGD rate mom (reg * 10)
  sgdUpdateLearningParamters o                     = o

data MnistOpts = MnistOpts FilePath FilePath Int Bool (Optimizer 'SGD) (Optimizer 'Adam)

mnist' :: Parser MnistOpts
mnist' = MnistOpts <$> argument str (metavar "TRAIN")
                   <*> argument str (metavar "VALIDATE")
                   <*> option auto (long "iterations" <> short 'i' <> value 15)
                 <*> flag False True (long "use-adam" <> short 'a')
                 <*> (OptSGD
                       <$> option auto (long "train_rate" <> short 'r' <> value 0.01)
                       <*> option auto (long "momentum" <> value 0.9)
                       <*> option auto (long "l2" <> value 0.0005)
                       )
                 <*> (OptAdam
                       <$> option auto (long "alpha" <> short 'r' <> value 0.001)
                       <*> option auto (long "beta1" <> value 0.9)
                       <*> option auto (long "beta2" <> value 0.999)
                       <*> option auto (long "epsilon" <> value 1e-4)
                       <*> option auto (long "lambda" <> value 1e-3)
                      )

main :: IO ()
main = do
    MnistOpts mnist vali iter useAdam sgd adam <- execParser (info (mnist' <**> helper) idm)
    putStrLn "Training convolutional neural network..."
    res <- if useAdam
      then runExceptT $ convTest iter mnist vali adam
      else runExceptT $ convTest iter mnist vali sgd

    case res of
      Right () -> pure ()
      Left err -> putStrLn err

readMNIST :: FilePath -> ExceptT String IO [(S ('D2 28 28), S ('D1 10))]
readMNIST mnist = ExceptT $ do
  mnistdata <- T.readFile mnist
  return $ traverse (A.parseOnly parseMNIST) (tail $ T.lines mnistdata)

parseMNIST :: A.Parser (S ('D2 28 28), S ('D1 10))
parseMNIST = do
  Just lab <- oneHot <$> A.decimal
  pixels   <- many (A.char ',' >> A.double)
  image    <- maybe (fail "Parsed row was of an incorrect size") pure (fromStorable . V.fromList $ map realToFrac pixels)
  return (image, lab)

max_stars_repo_name: schnecki/grenade

max_issues_repo_name: schnecki/grenade

max_forks_repo_name: schnecki/grenade


text:
import Control.Exception (evaluate)
import Data.Complex
import Test.Hspec
import Test.QuickCheck

import CommonFunctions
import JuliaSet

main :: IO ()
main = hspec $ do
  -- test CommonFunctions
  describe "CommonFunctions.f1" $ do
    it "should work with only real numbers" $
      f1 1 2 `shouldBe` 5

    it "should work with complex constant and real variable" $
      f1 (1 :+ 1) 2 `shouldBe` (5 :+ 1)

    it "should work with real constant and complex variable" $
      f1 2 (2 :+ 3) `shouldBe` ((-3) :+ 12)

    it "should work with only complex numbers" $
      f1 (7 :+ 9.25) (3 :+ 2) `shouldBe` (12 :+ 21.25)

    it "should use partial application to create function with constant c" $ do
      let f2 = f1 (16 :+ 4)
      f2 5 `shouldBe` (41 :+ 4)
      f2 5.5 `shouldBe` (46.25 :+ 4)
      f2 6 `shouldBe` (52 :+ 4)

  -- test JuliaSet
  describe "JuliaSet.julia" $ do
    let f = f1 ((-0.4) :+ 0.65)
    it "should produce 1 for z = -1.5 + 1.0i where f = z^2 -0.4 + 0.65i" $
      julia f ((-1.5) :+ (1.0)) 2 1000 `shouldBe` 1
    it "should be > 1 for z = -0.4 + 0.4i where f = z^2 -0.4 + 0.65i" $
      julia f ((-0.4) :+ 0.4) 2 1000 `shouldSatisfy` (> 1)


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

module Main where

import Control.Exception
import qualified Data.Array.Accelerate as A
import qualified Data.Array.Accelerate.Interpreter as ALI
import Data.Array.Accelerate.Math.Hilbert
import Data.Array.Accelerate.Math.Qtfd
import Data.Array.Accelerate.Math.WindowFunc
import Data.Attoparsec.Text
import ParseArgs
import ParseFile
import System.Directory
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
import qualified Data.Array.Accelerate.LLVM.Native as ALN
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
import qualified Data.Array.Accelerate.LLVM.PTX as ALP
#endif
import Data.Array.Accelerate.Array.Sugar as S
import qualified Data.ByteString.Builder as BB
import Data.Complex
import qualified Data.Double.Conversion.Convertable as DC
import Data.List
import Data.Monoid
import qualified Data.Text as T
import qualified Data.Text.IO as TI
import qualified Data.Vector as V
import GHC.Float
import Math.Gamma
import System.IO

main :: IO ()
main = do
  args <- opts
  makeQTFDAll args

-- | Transform data to text(using Float conversion. It is much faster) and write it to file.

writeData ::
  Handle            -- ^ output file
  -> [Float]       -- ^ Data to write
  -> Int            -- ^ number of element
  -> Int            -- ^ elements per row
  -> IO ()
writeData _ [] _ _  = return ()
writeData file (x:xs) w n = do
  TI.hPutStr file $ DC.toPrecision 5 x <> " "
  if (n `mod` w) == 0
  then do
    TI.hPutStr file "\n"
    writeData file xs w (n + 1)
  else writeData file xs w (n + 1)

makeString ::
     Handle
  -> [Float]       -- ^ Data to write
  -> Int            -- ^ number of element
  -> Int            -- ^ elements per row
  -> BB.Builder
  -> IO ()
makeString _ [] _ _ _  = return ()
makeString file (x:xs) w n acc =
  if (n `mod` w) == 0
  then BB.hPutBuilder file (acc <> DC.toPrecision 5 x <> " " <> "\n") >> makeString file xs w (n+1) ""
  else makeString file xs w (n+1) (acc <> DC.toPrecision 5 x <> " ")

-- | make Pseudo Wigner-Ville transform to all files in a directory.

makeQTFDAll :: Opts -> IO ()
makeQTFDAll opts = do
  let path = getPath opts
      nosAVGflag = subAVGflag opts
  files <- listDirectory path
  setCurrentDirectory path
  let fFiles = filter (isSuffixOf ".txt") files
      sAVGflag = A.constant $ not nosAVGflag
  case opts of
    (OptsPWV _ dev wlen wfun _) ->
      do
        if (even wlen)
        then error "Window length maust be odd !"
        else
          do
            let window = makeWindow wfun (A.unit $ A.constant wlen) :: A.Acc (A.Array A.DIM1 Float)
            mapM_ (makePWV dev window sAVGflag) fFiles
    (OptsWV path dev nosAVGflag) -> mapM_ (makeWV dev sAVGflag) fFiles
    (OptsCW path dev sigma uWindow uWindowFunc nWindow nWindowFunc normalise nosAVGflag) ->
      do
        if (even uWindow || even nWindow)
        then error "Window length maust be odd !"
        else
          let uWindowArray = makeWindow uWindowFunc (A.unit $ A.constant uWindow) :: A.Acc (A.Array A.DIM1 Float)
              nWindowArray = makeWindow nWindowFunc (A.unit $ A.constant nWindow) :: A.Acc (A.Array A.DIM1 Float)
              bothRectangular = uWindowFunc == Rect && nWindowFunc == Rect
          in mapM_ (makeCW dev sAVGflag sigma uWindowArray nWindowArray normalise bothRectangular) fFiles
    (OptsBJ path dev alpha uWindow uWindowFunc nWindow nWindowFunc nosAVGflag) ->
      do
        if (even uWindow || even nWindow)
        then error "Window length maust be odd !"
        else
          let uWindowArray = makeWindow uWindowFunc (A.unit $ A.constant uWindow) :: A.Acc (A.Array A.DIM1 Float)
              nWindowArray = makeWindow nWindowFunc (A.unit $ A.constant nWindow) :: A.Acc (A.Array A.DIM1 Float)
              bothRectangular = uWindowFunc == Rect && nWindowFunc == Rect
          in mapM_ (makeBJ dev sAVGflag alpha uWindowArray nWindowArray bothRectangular) fFiles
    (OptsEMBD path dev alpha beta uWindow uWindowFunc nWindow nWindowFunc normalise nosAVGflag) ->
      do
        if (even uWindow || even nWindow)
        then error "Window length maust be odd !"
        else let gammas = makeGammaArray beta uWindow
                 uWindowArray = makeWindow uWindowFunc (A.unit $ A.constant uWindow) :: A.Acc (A.Array A.DIM1 Float)
                 nWindowArray = makeWindow nWindowFunc (A.unit $ A.constant nWindow) :: A.Acc (A.Array A.DIM1 Float)
                 bothRectangular = uWindowFunc == Rect && nWindowFunc == Rect
             in mapM_ (makeEMBD dev sAVGflag alpha gammas uWindowArray nWindowArray bothRectangular normalise) fFiles
  return ()

makeCW :: CalcDev                   -- ^ Calculation device - CPU or GPU
  -> A.Exp Bool                     -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float                          -- ^ sigma
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if the core array should be normalised
  -> Bool                           -- ^ if Both windows are rectangular
  -> FilePath                       -- ^ Name of file with data.
  -> IO ()
makeCW dev sAVGflag sigma uWindowArray nWindowArray normalise bothRect file = do
  putStrLn $ "processing " ++ file ++ "..."
  text <- TI.readFile file
  let parseRes = parseOnly parseFile text
  case parseRes of
    (Left errStr) -> error $ errStr
    (Right dataF) -> mapM_ (startCW dev file sAVGflag sigma uWindowArray nWindowArray normalise bothRect) dataF

makeBJ :: CalcDev                   -- ^ Calculation device - CPU or GPU
  -> A.Exp Bool                     -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float                          -- ^ sigma
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if Both windows are rectangular
  -> FilePath                       -- ^ Name of file with data.
  -> IO ()
makeBJ dev sAVGflag alpha uWindowArray nWindowArray bothRect file = do
  putStrLn $ "processing " ++ file ++ "..."
  text <- TI.readFile file
  let parseRes = parseOnly parseFile text
  case parseRes of
    (Left errStr) -> error $ errStr
    (Right dataF) -> mapM_ (startBJ dev file sAVGflag alpha uWindowArray nWindowArray bothRect) dataF

-- | Make Pseudo Wigner-Ville transform to all columns in the given file

makePWV ::
  CalcDev                            -- ^ Calculation device - CPU or GPU
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> A.Exp Bool                      -- ^ Apply subtraction of the mean value from all elements in a column.
  -> FilePath                        -- ^ Name of file with data.
  -> IO ()
makePWV dev window sAVGflag file = do
  putStrLn $ "processing " ++ file ++ "..."
  text <- TI.readFile file
  let parseRes = parseOnly parseFile text
  case parseRes of
    (Left errStr) -> error $ errStr
    (Right dataF) -> mapM_ (startPWV dev window file sAVGflag) dataF

--makeEMBD dev alpha gammas uWindowArray nWindowArray bothRectangular normalise

makeEMBD :: CalcDev                   -- ^ Calculation device - CPU or GPU
  -> A.Exp Bool                     -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float                          -- ^ alpha
  -> [Float]                        -- ^ gammas
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if the core array should be normalised
  -> Bool                           -- ^ if Both windows are rectangular
  -> FilePath                       -- ^ Name of file with data.
  -> IO ()
makeEMBD dev sAVGflag alpha gammas uWindowArray nWindowArray normalise bothRect file = do
  putStrLn $ "processing " ++ file ++ "..."
  text <- TI.readFile file
  let parseRes = parseOnly parseFile text
  case parseRes of
    (Left errStr) -> error $ errStr
    (Right dataF) -> mapM_ (startEMBD dev file sAVGflag alpha gammas uWindowArray nWindowArray normalise bothRect) dataF

-- | Make Pseudo Wigner-Ville transform to the given column.
-- At first it makes subtraction of average value (if supAVGflag) and applies hilbert transform to make analytic signal. After that it applies Pseudo-Wigner transftorm
-- and save data to file with

startPWV ::
  CalcDev                          -- ^ Calculation device - CPU or GPU
  -> A.Acc (A.Array A.DIM1 Float) -- ^ Smoothing window in time-domain. Length must be odd.
  -> FilePath                      -- ^ Name of file with data.
  -> A.Exp Bool                    -- ^ Apply subtraction of the mean value from all elements in a column.
  -> (T.Text,[Float])             -- ^ Name of column and parsed data from column
  -> IO ()
startPWV dev window oldName sAVGflag (column_name,dataF) = do
  let newFName = oldName ++ "-" ++ T.unpack column_name ++ ".txt"
  putStr $ "   Creating file " ++ newFName ++ " ..."
  let leng = length dataF
      pData = A.fromList (A.Z A.:. leng) dataF
      appPWV = (pWignerVille window) . hilbert . supAVG sAVGflag
      processed = case dev of
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
                    CPU -> ALN.run1 appPWV pData
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
                    GPU -> ALP.run1 appPWV pData
#endif
                    CPU -> ALI.run1 appPWV pData
                    GPU -> error "Compiled without GPU support"
      pList = S.toList $! processed
  file <- openFile newFName WriteMode
  onException (writeData file pList leng 1) (removeFile newFName)
  hClose file
  putStrLn "Done !"

-- | Make Wigner-Ville transform to all columns in the given file

makeWV ::
  CalcDev        -- ^ Calculation device - CPU or GPU
  -> A.Exp Bool  -- ^ Apply subtraction of the mean value from all elements in a column.
  -> FilePath    -- ^ Name of file with data.
  -> IO ()
makeWV dev sAVGflag file = do
  putStrLn $ "processing " ++ file ++ "..."
  text <- TI.readFile file
  let parseRes = parseOnly parseFile text
  case parseRes of
    (Left errStr) -> error $ errStr
    (Right dataF) -> mapM_ (startWV dev file sAVGflag) dataF

-- | Make Wigner-Ville transform to the given column.
-- At first it makes subtraction of average value (if supAVGflag) and applies hilbert transform to make analytic signal. After that it applies Pseudo-Wigner transftorm
-- and save data to file with

startWV ::
  CalcDev                 -- ^ Calculation device - CPU or GPU
  -> FilePath             -- ^ Name of file with data.
  -> A.Exp Bool           -- ^ Apply subtraction of the mean value from all elements in a column.
  -> (T.Text,[Float])    -- ^ Name of column and parsed data from column
  -> IO ()
startWV dev oldName sAVGflag (column_name,dataF) = do
  let newFName = oldName ++ "-" ++ T.unpack column_name ++ ".txt"
  putStr $ "   Creating file " ++ newFName ++ " ..."
  let leng = length dataF
      pData = A.fromList (A.Z A.:. leng) dataF
      appWV = wignerVilleNew . hilbert . supAVG sAVGflag
      appWVGPU = wignerVille . hilbert . supAVG sAVGflag
      processed = case dev of
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
                    CPU -> ALN.run1 appWV pData
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
                    GPU -> ALP.run1 appWV pData
#endif
                    CPU -> ALI.run1 appWV pData
                    GPU -> error "Compiled without GPU support"
      pList = S.toList $  processed
  file <- openFile newFName WriteMode
  hSetBinaryMode file True >> hSetBuffering file LineBuffering
  -- let resultBsBuilder = makeString pList leng 1 ""
  onException (makeString file pList leng 1 "") (removeFile newFName) >> hFlush file
  hClose file
  putStrLn "Done !"

startCW ::
  CalcDev                 -- ^ Calculation device - CPU or GPU
  -> FilePath             -- ^ Name of file with data.
  -> A.Exp Bool           -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float               -- ^ Sigma
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if the core array should be normalised
  -> Bool                           -- ^ if Both windows are rectangular
  -> (T.Text,[Float])    -- ^ Name of column and parsed data from column
  -> IO ()
startCW dev oldName sAVGflag sigma uWindowArray nWindowArray normalise bothRect (column_name,dataF) = do
  let newFName = oldName ++ "-" ++ T.unpack column_name ++ ".txt"
  putStr $ "   Creating file " ++ newFName ++ " ..."
  let leng = length dataF
      uWindow = A.length uWindowArray
      nWindow = A.length nWindowArray
      mWindowArrays = if bothRect then Nothing else Just (uWindowArray,nWindowArray)
      pData = A.fromList (A.Z A.:. leng) dataF
      appCW = (choiWilliams (A.constant sigma) mWindowArrays uWindow nWindow normalise) . hilbert . supAVG sAVGflag
      processed = case dev of
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
                    CPU -> ALN.run1 appCW pData
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
                    GPU -> ALP.run1 appCW pData
#endif
                    CPU -> ALI.run1 appCW pData
                    GPU -> error "Compiled without GPU support"
      pList = S.toList $  processed
  file <- openFile newFName WriteMode
  onException (writeData file pList leng 1) (removeFile newFName)
  hClose file
  putStrLn "Done !"

startBJ ::
  CalcDev                 -- ^ Calculation device - CPU or GPU
  -> FilePath             -- ^ Name of file with data.
  -> A.Exp Bool           -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float               -- ^ Sigma
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if Both windows are rectangular
  -> (T.Text,[Float])    -- ^ Name of column and parsed data from column
  -> IO ()
startBJ dev oldName sAVGflag alpha uWindowArray nWindowArray bothRect (column_name,dataF) = do
  let newFName = oldName ++ "-" ++ T.unpack column_name ++ ".txt"
  putStr $ "   Creating file " ++ newFName ++ " ..."
  let leng = length dataF
      uWindow = A.length uWindowArray
      nWindow = A.length nWindowArray
      mWindowArrays = if bothRect then Nothing else Just (uWindowArray,nWindowArray)
      pData = A.fromList (A.Z A.:. leng) dataF
      appCW = (bornJordan (A.constant alpha) mWindowArrays uWindow nWindow) . hilbert . supAVG sAVGflag
      processed = case dev of
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
                    CPU -> ALN.run1 appCW pData
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
                    GPU -> ALP.run1 appCW pData
#endif
                    CPU -> ALI.run1 appCW pData
                    GPU -> error "Compiled without GPU support"
      pList = S.toList $  processed
  file <- openFile newFName WriteMode
  onException (writeData file pList leng 1) (removeFile newFName)
  hClose file
  putStrLn "Done !"

startEMBD ::
  CalcDev                 -- ^ Calculation device - CPU or GPU
  -> FilePath             -- ^ Name of file with data.
  -> A.Exp Bool           -- ^ Apply subtraction of the mean value from all elements in a column.
  -> Float                -- ^ Alpha
  -> [Float]              -- ^ Gammas
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in frequensy-domain. Length must be odd.
  -> A.Acc (A.Array A.DIM1 Float)   -- ^ Smoothing window in time-domain. Length must be odd.
  -> Bool                           -- ^ if the core array should be normalised
  -> Bool                           -- ^ if Both windows are rectangular
  -> (T.Text,[Float])    -- ^ Name of column and parsed data from column
  -> IO ()
startEMBD dev oldName sAVGflag alpha gammas uWindowArray nWindowArray normalise bothRect (column_name,dataF) = do
  let newFName = oldName ++ "-" ++ T.unpack column_name ++ ".txt"
  putStr $ "   Creating file " ++ newFName ++ " ..."
  let leng = length dataF
      uWindow = A.length uWindowArray
      nWindow = A.length nWindowArray
      gammasArr = A.use $ A.fromList (A.Z A.:. (length gammas)) gammas
      mWindowArrays = if bothRect then Nothing else Just (uWindowArray,nWindowArray)
      pData = A.fromList (A.Z A.:. leng) dataF
      appEMBD = (eModifiedB (A.constant alpha) gammasArr mWindowArrays uWindow nWindow normalise) . hilbert . supAVG sAVGflag
      processed = case dev of
#ifdef ACCELERATE_LLVM_NATIVE_BACKEND
                    CPU -> ALN.run1 appEMBD pData
#endif
#ifdef ACCELERATE_LLVM_PTX_BACKEND
                    GPU -> ALP.run1 appEMBD pData
#endif
                    CPU -> ALI.run1 appEMBD pData
                    GPU -> error "Compiled without GPU support"
      pList = S.toList $  processed
  file <- openFile newFName WriteMode
  onException (writeData file pList leng 1) (removeFile newFName)
  hClose file
  putStrLn "Done !"

-- | Subtraction of average value from all elements of column

supAVG :: A.Exp Bool -> A.Acc (Array DIM1 Float) -> A.Acc (Array DIM1 Float)
supAVG flag arr =
  A.acond flag (A.map (\x -> x - avg) arr) arr
  where leng = A.length arr
        avg = (A.the $ A.sum arr)/(A.fromIntegral leng)

makeGammaArray :: Float -> Int -> [Float]
makeGammaArray beta wLength =
  let dwLength = fromIntegral wLength :: Float
      v = [(-0.5), ((-0.5) + 1.0/dwLength)..0.5]
      s = map (\x -> beta :+ pi*x) v
      f = map (\x -> abs $ ((magnitude $ gamma x) ^ 2)/((gamma beta) ^ 2)) s
      n = floor $ (fromIntegral $ length f)/2.0
      frst = Data.List.take (n-1) f
      secnd  = Data.List.drop n f
  in secnd ++ frst

max_stars_repo_name: Haskell-mouse/WVille

max_issues_repo_name: Haskell-mouse/WVille

max_forks_repo_name: Haskell-mouse/WVille


text:
-- | An implementation of ResourceHandle-based variables.
--
-- The main difference between this and 'Ref'-based variables is
-- that reads are explicit, via the 'readValue' op.
--
-- TODO: given that distinction, figure out a good story around
-- gradients and save/restore.  Then, merge this module into
-- TensorFlow.Ops.
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE RecursiveDo #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE OverloadedStrings #-}
module TensorFlow.Variable
    ( Variable
    , variable
    , variable'
    , readValue
    , initializedValue
    , initializedVariable
    , initializedVariable'
    , zeroInitializedVariable
    , zeroInitializedVariable'
    , assign
    , assign'
    , assignAdd
    , assignAdd'
    , resourceApplyAdam
    , resourceApplyAdam'
    ) where

import qualified Data.Complex
import qualified Data.Int
import qualified Data.Word
import Data.Text.Encoding (encodeUtf8)
import Lens.Family2 ((.~), (&))
import TensorFlow.Core
import TensorFlow.Build (opDef)
import TensorFlow.BuildOp (buildInputs, pureOp, OpParams)
import TensorFlow.Output (opInputs, unNodeName)
import TensorFlow.Tensor (Rendered(..), ToTensor(..), renderValue, tensorNodeName)
import TensorFlow.Types (tensorType)
import qualified TensorFlow.GenOps.Core as CoreOps
import TensorFlow.Ops (zeros)

data Variable a = Variable
    { variableHandle   :: Tensor Value ResourceHandle
    , initializedValue :: Maybe (Tensor Value a)
      -- ^ The initial value of a 'Variable' created with 'initializedVariable'.
    }

instance Rendered Variable where
    renderedOutput = renderedOutput . variableHandle

instance ToTensor Variable where
    toTensor = readValue

-- | Creates a new, uninitialized variable.
variable :: (MonadBuild m, TensorType a) => Shape -> m (Variable a)
variable = variable' id

variable' :: forall m a . (MonadBuild m, TensorType a)
                    => OpParams -> Shape -> m (Variable a)
variable' params s = build $ do
    -- Each variable needs a unique "shared_name".  Use MonadFix to
    -- set the attribute to the same name as the variable itself, without
    -- exposing more internals of the Build module.
    rec t <- CoreOps.varHandleOp' (params . (opAttr "shared_name" .~ n))
                                    (tensorType (undefined :: a)) s
        let n = encodeUtf8 $ unNodeName $ tensorNodeName t
    return $ Variable t Nothing

-- | Creates a variable initialized to the given value.
-- Initialization happens next time session runs.
initializedVariable :: (MonadBuild m, TensorType a)
                    => Tensor v a -> m (Variable a)
initializedVariable = initializedVariable' id

initializedVariable' :: forall a m v . (MonadBuild m, TensorType a)
                    => OpParams -> Tensor v a -> m (Variable a)
initializedVariable' params initializer = do
    -- The shape is not known initially.
    (Variable h Nothing :: Variable a) <- variable' params (Shape [])
    initializer' <- renderValue initializer
    i <- CoreOps.assignVariableOp h initializer'
    addInitializer =<< group i
    return (Variable h (Just initializer'))

-- | Creates a zero-initialized variable with the given shape.
zeroInitializedVariable
  :: (MonadBuild m, TensorType a, Num a) => Shape -> m (Variable a)
zeroInitializedVariable = zeroInitializedVariable' id

zeroInitializedVariable'
  :: (MonadBuild m, TensorType a, Num a) => OpParams -> Shape -> m (Variable a)
zeroInitializedVariable' params = initializedVariable' params . zeros

-- | Gets the value stored in a variable.
--
-- Note that this op is stateful since it depends on the value of the variable;
-- however, it may be CSE'd with other reads in the same context.  The context can
-- be fixed by using 'render' along with (for example) 'withControlDependencies'.
-- For example:
--
-- >   runSession $ do
-- >     v <- variable []
-- >     a <- assign v 24
-- >     r <- withControlDependencies a $ render $ readValue v + 18
-- >     result <- run r
-- >     liftIO $ (42 :: Float) @=? unScalar result
--
--
readValue :: TensorType a => Variable a -> Tensor Build a
readValue = readValue' id

readValue' :: forall a . TensorType a
    => OpParams -> Variable a -> Tensor Build a
readValue' params (Variable h _)
    = pureOp [] $ do
        os <- buildInputs h
        pure $ opDef "ReadVariableOp"
                & (params
                    . (opAttr "dtype" .~ tensorType (undefined :: a))
                    . (opInputs .~ os))

-- | Sets the value of a variable.
assign :: (MonadBuild m, TensorType a)
    => Variable a -> Tensor v a -> m ControlNode
assign = assign' id

assign' :: (MonadBuild m, TensorType a)
    => OpParams -> Variable a -> Tensor v a -> m ControlNode
assign' params (Variable h _) v = CoreOps.assignVariableOp' params h v

-- | Increments the value of a variable.
assignAdd :: (MonadBuild m, TensorType a)
    => Variable a -> Tensor v a -> m ControlNode
assignAdd = assignAdd' id

assignAdd' :: (MonadBuild m, TensorType a)
    => OpParams -> Variable a -> Tensor v a -> m ControlNode
assignAdd' params (Variable h _) v = CoreOps.assignAddVariableOp' params h v

-- | Update '*var' according to the Adam algorithm.
--
-- lr_t <- learning_rate * sqrt(1 - beta2^t) / (1 - beta1^t)
-- m_t <- beta1 * m_{t-1} + (1 - beta1) * g_t
-- v_t <- beta2 * v_{t-1} + (1 - beta2) * g_t * g_t
-- variable <- variable - lr_t * m_t / (sqrt(v_t) + epsilon)
resourceApplyAdam ::
    (MonadBuild m,
     OneOf '[(Data.Complex.Complex Double),
             (Data.Complex.Complex Float),
             Data.Int.Int16,
             Data.Int.Int32,
             Data.Int.Int64, Data.Int.Int8,
             Data.Word.Word16,
             Data.Word.Word8, Double,
             Float] t)
    => Variable t -- ^ __var__: Should be from a Variable().
    -> Variable t -- ^ __m__: Should be from a Variable().
    -> Variable t -- ^ __v__: Should be from a Variable().
    -> Tensor v1 t -- ^ __beta1_power__: Must be a scalar.
    -> Tensor v2 t -- ^ __beta2_power__: Must be a scalar.
    -> Tensor v3 t -- ^ __lr__: Scaling factor. Must be a scalar.
    -> Tensor v4 t -- ^ __beta1__: Momentum factor. Must be a scalar.
    -> Tensor v5 t -- ^ __beta2__: Momentum factor. Must be a scalar.
    -> Tensor v6 t -- ^ __epsilon__: Ridge term. Must be a scalar.
    -> Tensor v7 t -- ^ __grad__: The gradient.
    -> m (ControlNode)
resourceApplyAdam = resourceApplyAdam' id

resourceApplyAdam' ::
    (MonadBuild m,
     OneOf '[(Data.Complex.Complex Double),
             (Data.Complex.Complex Float),
             Data.Int.Int16, Data.Int.Int32,
             Data.Int.Int64, Data.Int.Int8,
             Data.Word.Word16, Data.Word.Word8, Double,
             Float] t)
    => OpParams
    -> Variable t -- ^ __var__: Should be from a Variable().
    -> Variable t -- ^ __m__: Should be from a Variable().
    -> Variable t -- ^ __v__: Should be from a Variable().
    -> Tensor v1 t -- ^ __beta1_power__: Must be a scalar.
    -> Tensor v2 t -- ^ __beta2_power__: Must be a scalar.
    -> Tensor v3 t -- ^ __lr__: Scaling factor. Must be a scalar.
    -> Tensor v4 t -- ^ __beta1__: Momentum factor. Must be a scalar.
    -> Tensor v5 t -- ^ __beta2__: Momentum factor. Must be a scalar.
    -> Tensor v6 t -- ^ __epsilon__: Ridge term. Must be a scalar.
    -> Tensor v7 t -- ^ __grad__: The gradient.
    -> m (ControlNode)
resourceApplyAdam' params (Variable var _) (Variable m _) (Variable v _) =
    CoreOps.resourceApplyAdam' params var m v

max_stars_repo_name: Qinka/tensor-haskell

max_issues_repo_name: Qinka/tensor-haskell

max_forks_repo_name: Qinka/tensor-haskell


text:
import Test.HUnit
import FFT.FFTControl
import FFT.Orig
import FFT.Samples
import System.Environment
import Data.Complex

testFFTbflySPipeline = TestCase (assertEqual "FFT-bflyS-Pipline" expected actual)
    where
        expected = (sum $ fftbflySPipeline samples512)
        actual = (sum $ fft samples512)
        samples512 = samples 1 512

testFFTDc = TestCase (assertEqual "DFTDc" expected actual)
    where
        expected = (sum $ fftdc samples512)
        actual = (sum $ fft samples512)
        samples512 = samples 1 512

tests = TestList [TestLabel "testFFTbFlysPipeline" testFFTbflySPipeline,
		TestLabel "testFFTDc" testFFTDc
                    ]

main = runTestTT tests

max_stars_repo_name: BNJHope/parallel-fourier-transforms

max_issues_repo_name: BNJHope/parallel-fourier-transforms

max_forks_repo_name: BNJHope/parallel-fourier-transforms


text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE TemplateHaskell   #-}
{-# LANGUAGE CApiFFI           #-}
-- |
-- Copyright: 2018 (C) Tweag I/O Limited.
--
-- inline-c context.
module Foreign.R.Context
  ( rCtx
  , SEXPREC
  , SEXP0
  , Logical(..)
  ) where

import Data.Complex
import qualified Data.Map as Map
import Foreign.C
import Foreign.Ptr
import Foreign.Storable
import Language.C.Types (TypeSpecifier(TypeName))
import Language.C.Inline.Context (Context(..))
import Internal.Error

#include <Rinternals.h>

data {-# CTYPE  "SEXPREC" #-} SEXPREC

type SEXP0 = Ptr SEXPREC

-- | R uses three-valued logic.
data {-# CTYPE "Logical" #-} Logical = FALSE
             | TRUE
             | NA
-- XXX no Enum instance because NA = INT_MIN, not representable as an Int on
-- 32-bit systems.
               deriving (Eq, Ord, Show)

instance Storable Logical where
  sizeOf _       = sizeOf (undefined :: CInt)
  alignment _    = alignment (undefined :: CInt)
  poke ptr FALSE = poke (castPtr ptr) (0 :: CInt)
  poke ptr TRUE  = poke (castPtr ptr) (1 :: CInt)
  -- Currently NA_LOGICAL = INT_MIN.
  poke ptr NA    = poke (castPtr ptr) (#{const INT_MIN} :: CInt)
  peek ptr = do
      x <- peek (castPtr ptr)
      case x :: CInt of
          0 -> return FALSE
          1 -> return TRUE
          #{const INT_MIN} -> return NA
          _ -> failure "Storable Logical peek" "Not a Logical."

rCtx :: Context
rCtx = mempty { ctxTypesTable = Map.fromList tytabs }
  where
    tytabs =
      [ (TypeName "SEXP", [t| Ptr SEXPREC |])
      , (TypeName "Rcomplex", [t| Complex Double |])
      ]

ext: hsc

max_stars_repo_name: reactormonk/HaskellR

max_issues_repo_name: reactormonk/HaskellR

max_forks_repo_name: reactormonk/HaskellR


text:
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE FlexibleContexts #-}

module Data.AdaptiveGeneticAlgorithm
  ( adaptiveGeneticAlgorithm )
  where

import Control.Monad.State.Strict
import Data.AdaptiveCoordinateDescent.Internal
import Data.List ( sortBy )
import Data.Ord
import Data.Foldable
import Data.Traversable
import qualified Data.Vector as V
import Numeric.LinearAlgebra hiding ( toList )
import Pipes hiding ( for )
import qualified Pipes.Prelude as P
import System.Random.MWC

data GeneticAlgorithmState f = GeneticAlgorithmState
  { principalComponents :: (Vector Double, [Vector Double])
  , populationPool :: [f Double]
  , rng :: GenIO }

generateCauchy :: GenIO -> Double -> Double -> IO Double
generateCauchy rng mean scale = do
  r <- uniformR (0, 1) rng
  return $ (tan (pi * (r - 0.5)) * scale) + mean

principalComponentsFromPopulation :: Foldable f => [f Double] -> (Vector Double, [Vector Double])
principalComponentsFromPopulation population =
  let mat = normalize $ (length population><num_coordinates)
                        (mconcat $ fmap toList population)
      (evec, mat2) = rightSV mat

   in (evec, toColumns mat2)
 where
  num_coordinates = length $ toList $ head population

sortByFitness :: Monad m
              => (f Double -> m Double)
              -> [f Double]
              -> m [(f Double, Double)]
sortByFitness evaluate items = do
  evaluated_items <- for items $ \item -> (,) item <$> evaluate item
  return $ sortBy (comparing snd) evaluated_items

adaptiveGeneticAlgorithm :: (MonadIO m, Traversable f)
                         => (f Double -> m Double) -- ^ Cost function. Smaller cost = better model.
                         -> f Double
                         -> Int
                         -> Producer (f Double, Double) m void
adaptiveGeneticAlgorithm evaluate initial_params how_many_to_keep = do
  rrng <- liftIO createSystemRandom

  initial_items <- liftIO $ P.toListM (makePerturbed initial_params rrng (how_many_to_keep*2))

  best_initial_items' <- lift $ fmap (take how_many_to_keep) $ sortByFitness evaluate initial_items

  yield (head best_initial_items')

  let best_initial_items = fmap fst best_initial_items'
      pcomponents = principalComponentsFromPopulation best_initial_items

  evalStateT piper (GeneticAlgorithmState { principalComponents = pcomponents
                                          , rng = rrng
                                          , populationPool = best_initial_items })
 where
  piper = do
    current_pool <- V.fromList . populationPool <$> get
    (pcomponent_scale_vec, pcomponents) <- principalComponents <$> get
    for_ (zip [0..] pcomponents) $ \(index, pcomponent) -> do
      let pcomponent_scale = pcomponent_scale_vec `atIndex` index
      parent_index <- randomI (0, V.length current_pool-1)
      let parent = current_pool V.! parent_index
      scalef <- randomCauchy
      let final_scale = scalef * pcomponent_scale
      addCandidate (plus parent (fromVector (final_scale `scale` pcomponent) parent))

    current_candidates <- populationPool <$> get
    best_candidates <- lift $ lift $ sortByFitness evaluate current_candidates

    let truncated_candidates = take how_many_to_keep best_candidates
        pcomponents = principalComponentsFromPopulation $ fmap fst truncated_candidates

    modify $ \old -> old { principalComponents = pcomponents
                         , populationPool = fmap fst truncated_candidates }

    lift $ yield (head best_candidates)

    piper

  randomCauchy = do
    rrng <- rng <$> get
    liftIO $ generateCauchy rrng 0 1

  randomI range = do
    rrng <- rng <$> get
    liftIO $ uniformR range rrng

  addCandidate item = do
    st <- get
    put $ st { populationPool = item:populationPool st }

max_stars_repo_name: Noeda/adaptivecoordinatedescent

max_issues_repo_name: Noeda/adaptivecoordinatedescent

max_forks_repo_name: Noeda/adaptivecoordinatedescent