text:
{-# LANGUAGE AllowAmbiguousTypes, RankNTypes, UnicodeSyntax, NoMonomorphismRestriction #-}
{-# LANGUAGE BangPatterns, FlexibleInstances #-}

module Optimizer where

import Utils
import Style
import GenOptProblem
import Numeric.AD
import Numeric.AD.Internal.On
import Numeric.AD.Internal.Reverse
import Numeric.AD.Internal.Sparse
import qualified Numeric.LinearAlgebra as L
import Debug.Trace
import System.Random
import System.Console.ANSI
import Data.List (foldl')

default (Int, Float)

------ Opt types, util functions, and params

type ObjFn1 a = forall a . (Autofloat a) => [a] -> a

-- used for duf
type ObjFn2 a = forall a . (Autofloat a) => [a] -> [a] -> a

type GradFn a = forall a . (Autofloat a) => [a] -> [a]

----- Various consts

nanSub :: (Autofloat a) => a
nanSub = 0

infinity :: Floating a => a
infinity = 1/0 -- x/0 == Infinity for any x > 0 (x = 0 -> Nan, x < 0 -> -Infinity)
-- all numbers are smaller than infinity except infinity, to which it's equal

----- Hyperparameters

weightGrowthFactor :: (Autofloat a) => a -- for EP weight
weightGrowthFactor = 10

epsUnconstr :: Floating a => a
epsUnconstr = 10 ** (-2)

epStop :: Floating a => a -- for EP diff
epStop = 10 ** (-5)
-- epStop = 60 ** (-3)
-- epStop = 10 ** (-1)
-- epStop = 0.05

-- Parameters for Armijo-Wolfe line search
-- NOTE: must maintain 0 < c1 < c2 < 1
c1 :: Floating a => a
c1 = 0.001 -- for Armijo, corresponds to alpha in backtracking line search (see below for explanation)
-- smaller c1 = shallower slope = less of a decrease in fn value needed = easier to satisfy
-- turn Armijo off: c1 = 0
-- Nocedal p38: "In practice, c1 is chosen to be quite small, say 10−4."

c2 :: Floating a => a
c2 = 0.9 -- for Wolfe, is the factor decrease needed in derivative value
-- new directional derivative value / old DD value <= c2
-- smaller c2 = smaller new derivative value = harder to satisfy
-- turn Wolfe off: c1 = 1 (basically backatracking line search only)
-- Nocedal p39: "Typical values of c2 are 0.9 when the search direction pk is chosen by a Newton or quasi-Newton method"

-- true = force linesearch halt if interval gets too small; false = no forced halt
intervalMin = True

useLineSearch :: Bool
useLineSearch = True

useAutodiff :: Bool
useAutodiff = True

constT :: Floating a => a
constT = 0.001

debugOpt = True
debugLineSearch = False
debugBfgs = True

trb :: String -> a -> a
trb s x = if debugBfgs then trace "---" $ trace s x else x -- prints in left to right order

tro :: String -> a -> a
tro s x = if debugOpt then trace "---" $ trace s x else x -- prints in left to right order

trl :: Show a => String -> a -> a
trl s x = if debugLineSearch then trace "---" $ trace s $ traceShowId x else x -- prints in left to right order

----- Convergence criteria

-- convergence criterion for EP
-- if you want to use it for UO, needs a different epsilon
epStopCond :: (Autofloat a) => [a] -> [a] -> a -> a -> Bool
epStopCond x x' fx fx' =
           tro ("EP: \n||x' - x||: " ++ (show $ norm (x -. x'))
           ++ "\n|f(x') - f(x)|: " ++ (show $ abs (fx - fx'))) $
           (norm (x -. x') <= epStop) || (abs (fx - fx') <= epStop)

unconstrainedStopCond :: (Autofloat a) => [a] -> Bool
unconstrainedStopCond gradEval = norm gradEval < epsUnconstr

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

-- Policies

stepPolicy :: State -> (Params, PolicyParams)
stepPolicy s =
    -- Check overall convergence first
    let epStatus = optStatus $ (paramsr s) in
    let pparams = policyParams s in
    case epStatus of

    -- Generate new objective function and replace the optimization and policy params accordingly
    EPConverged ->
                -- TODO: clean up the step incrementing
        let pparams' = pparams { policySteps = 1 + policySteps pparams } in
        let (policyRes, psNew) = (policyFn s) (objFns s) (constrFns s) pparams' in -- See what the policy function wants
        case policyRes of
            Nothing     -> (paramsr s, pparams' { policyState = psNew }) -- steps incremented, policy done

            Just newFns -> -- Policy keeps going
                let objFnNew = genObjfn (transr s) (filter isObjFn newFns) (filter isConstr newFns) (varyingPaths s)
                    -- TODO: check that these inputs are right
                    -- Change obj function and restart optimization
                    pparamsNew = pparams' { policyState = psNew,
                                            currFns = newFns }
                    paramsNew = Params { weight = initWeight,
                                         optStatus = NewIter,
                                         overallObjFn = objFnNew,
                                         bfgsInfo = defaultBfgsParams }
                in tro ("Step policy, EP converged, new params:\n" ++ show (paramsNew, pparamsNew, newFns)) $ (paramsNew, pparamsNew)

    -- If not converged, optimize as usual, don't change policy mid-optimization
    _ -> tro ("Step policy, EP not converged, new params:\n" ++ show (paramsr s, pparams)) $ (paramsr s, pparams)

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

-- Main optimization functions

step :: State -> State
step s = let (state', params') = stepShapes (oConfig s) (paramsr s) (varyingState s) (rng s)
             s'                = s { varyingState = state',
                                     paramsr = params' }
             -- NOTE: we intentionally discard the random generator here because
             -- we want to have consistent computation output in a single
             -- optimization session
             -- For the same reason, all subsequent step* functions such as
             -- stepShapes do not return the new random generator
             (!shapes', _, _)     = evalTranslation s'

             -- Check the state and see if the overall objective function should be changed
             -- The policy may change EPConverged to a new iteration before the frontend sees it
             (paramsNew, pparamsNew) = stepPolicy s'

             -- For debugging
             oldParams = paramsr s

         in tro ("Params: \n" ++ show oldParams ++ "\n:") $
            s' { shapesr = shapes',
                 paramsr = paramsNew,
                 policyParams = pparamsNew }
            -- note: trans is not updated in state

-- Note use of realToFrac to generalize type variables (on the weight and on the varying state)

-- implements exterior point algo as described on page 6 here:
-- https://www.me.utexas.edu/~jensen/ORMM/supplements/units/nlp_methods/const_opt.pdf
stepShapes :: OptConfig -> Params -> [Float] -> StdGen -> ([Float], Params)
stepShapes config params vstate g = -- varying state
         -- if null vstate then error "empty state in stepshapes" else
         let (epWeight, epStatus) = (weight params, optStatus params) in
         case epStatus of

         -- start the outer EP optimization and the inner unconstrained optimization, recording initial EPstate
         NewIter -> let status' = UnconstrainedRunning (map realToFrac vstate) in
                    (vstate', params { weight = initWeight, optStatus = status', bfgsInfo = defaultBfgsParams } )
         -- check *weak* convergence of inner unconstrained opt.
         -- if UO converged, set opt state to converged and update UO state (NOT EP state)
         -- if not, keep running UO (inner state implicitly stored)
         -- note convergence checks are only on the varying part of the state
         UnconstrainedRunning lastEPstate ->  -- doesn't use last EP state
           -- let unconstrConverged = unconstrainedStopCond gradEval in
           let unconstrConverged = epStopCond vstate vstate' (objFnApplied vstate) (objFnApplied vstate') in
               -- Two stopping conditions
               -- unconstrainedStopCond gradEval in
           if unconstrConverged then
              let status' = UnconstrainedConverged lastEPstate in -- update UO state only!
              (vstate', params { optStatus = status', bfgsInfo = defaultBfgsParams }) -- note vstate' (UO converged), not vstate
           else (vstate', params { bfgsInfo = bfgs' }) -- update UO state but not EP state; UO still running

         -- check EP convergence. if converged then stop, else increase weight, update states, and run UO again
         -- TODO some trickiness about whether unconstrained-converged has updated the correct state
         -- and whether to check WRT the updated state or not
         UnconstrainedConverged lastEPstate ->
           let epConverged = epStopCond lastEPstate (map r2f vstate) -- lastEPstate is last state for converged UO
                                   (objFnApplied lastEPstate) (objFnApplied (map r2f vstate)) in
           if epConverged then
              let status' = EPConverged in -- no more EP state
              (vstate, params { optStatus = status', bfgsInfo = defaultBfgsParams }) -- do not update UO state
           -- update EP state: to be the converged state from the most recent UO
           else let status' = UnconstrainedRunning (map realToFrac vstate) in -- increase weight
                let epWeight' = weightGrowthFactor * epWeight in
                -- trace ("Unconstrained converged. New weight: " ++ show epWeight') $
                      (vstate, params { weight = epWeight', optStatus = status', bfgsInfo = defaultBfgsParams })

         -- done; don't update obj state or params; user can now manipulate
         EPConverged -> (vstate, params { bfgsInfo = defaultBfgsParams } )

         -- TODO: implement EPConvergedOverride (for when the magnitude of the gradient is still large)

         -- TODO factor out--only unconstrainedRunning needs to run stepObjective, but EPconverged needs objfn
        where (vstate', gradEval, bfgs') = stepWithObjective config g params vstate
              objFnApplied = (overallObjFn params) g (r2f $ weight params)

-- Given the time, state, and evaluated gradient (or other search direction) at the point,
-- return the new state

stepT :: Float -> Float -> Float -> Float
stepT dt x dfdx = x - dt * dfdx

-- Calculates the new state by calculating the directional derivatives (via autodiff)
-- and timestep (via line search), then using them to step the current state.
-- Also partially applies the objective function.
stepWithObjective :: OptConfig -> StdGen -> Params -> [Float] -> ([Float], [Float], BfgsParams)
stepWithObjective config g params state =
          -- get timestep via line search, and evaluated gradient at the state
          let (t', gradEval, gradToUse, bfgs') = timeAndGrad config params objFnApplied state
              -- step each parameter of the state with the time and gradient
              state' = map (\(v, dfdv) -> stepT t' v dfdv) (zip state $ gradToUse)
              (fx, fx') = (objFnApplied state, objFnApplied state')
          in -- if fx' > fx then error ("Error: new energy is greater than old energy: " ++ show (fx', fx)) else
             tro ("\n----------------------------------------\n"
                   ++ "\nopt params: \n" ++ (show params)
                   ++ "\n||x' - x||: " ++ (show $ norm (state -. state'))
                   ++ "\n|f(x') - f(x)|: " ++
                  (show $ abs (fx' - fx))
                   ++ "\nf(x'): \n" ++ (show fx')
                   ++ "\ngradEval: \n" ++ (show gradEval)
                   ++ "\n||gradEval||: \n" ++ (show $ norm gradEval)
                   ++ "\ngradToUse: \n" ++ (show gradToUse)
                   ++ "\n||gradToUse||: \n" ++ (show $ norm gradToUse)
                   -- ++ "\nhessian: \n" ++ (show $ h)
                   ++ "\nbfgs': \n" ++ (show bfgs') -- TODO: use trb
                   ++ "\n timestep: \n" ++ (show t')
                   ++ "\n original state: \n" ++ (show state)
                   ++ "\n new state: \n" ++ (show state')
                  )
             (state', gradEval, bfgs')

          where objFnApplied :: ObjFn1 b
                objFnApplied = (overallObjFn params) g cWeight

                cWeight = r2f $ weight params
                -- realToFrac generalizes the type variable `a` to the type variable `b`, which timeAndGrad expects

-- a version of grad with a clearer type signature
appGrad :: (Autofloat a) => (forall b . (Autofloat b) => [b] -> b) -> [a] -> [a]
appGrad f l = grad f l

instance Show (Numeric.AD.Internal.On.On
                             (Numeric.AD.Internal.Reverse.Reverse
                                s (Numeric.AD.Internal.Sparse.Sparse a))) where
         show a = "error: not sure how to derive show for hessian element"

appHess :: (Autofloat a) => (forall b . (Autofloat b) => [b] -> b) -> [a] -> [[a]]
appHess f l = hessian f l

-- Precondition the gradient
gradP :: OptConfig -> BfgsParams -> [Double] -> ObjFn1 a -> [Float] -> ([Double], BfgsParams)
gradP config bfgsParams gradEval f state =
      let x_k = L.vector $ map r2f state
          grad_fx_k = L.vector gradEval
      in case optMethod config of
      GradientDescent -> (gradEval, bfgsParams)

      Newton -> -- Precondition gradient with the pseudoinverse of the hessian
          let h = appHess f state
              h_list = (map r2f $ concat h) :: [Double]
              hinv = L.pinv $ L.matrix (length gradEval) $ h_list
              gradPreconditioned = hinv L.#> (L.vector gradEval) in
          (L.toList gradPreconditioned, bfgsParams)

      BFGS -> -- Approximate inverse of hessian with the change in gradient (see Nocedal S9.1 p224)
           let grad_val = lastGrad bfgsParams
               h_val = invH bfgsParams
               state_val = lastState bfgsParams
           in case (h_val, grad_val, state_val) of
              (Nothing, Nothing, Nothing) -> -- First step. Initialize the approximation to the identity (not clear how else to approximate it)
                    -- k=0 steps from x_0 to x_1: so on the first step, we take a normal gradient descent step
                    let h_0 = L.ident $ length gradEval
                    in (gradEval, bfgsParams { lastState = Just x_k, lastGrad = Just grad_fx_k, invH = Just h_0 })

              (Just h_km1, Just grad_fx_km1, Just x_km1) -> -- For x_{k+1}, to compute H_k, we need the (k-1) info
                    -- Our convention is that we are always working "at" k to compute k+1
                    -- x_0 doesn't require any H; x_1 (the first step) with k = 0 requires H_0
                    -- x_2 (the NEXT step) with k=1 requires H_1. For example>
                    -- x_2 = x_1 - alpha_1 H_1 grad f(x_1)   [GD step]
                    -- H_1 = V_0 H_0 V_0 + rho_0 s_0 s_0^T   [This is confusing because the book adds an extra +1 to the H index]
                    -- V_0 = I - rho_0 y_0 s_0^T
                    -- rho_0 = 1 / y_0^T s_0
                    -- s_0 = x_1 - x_0
                    -- y_0 = grad f(x_1) - grad f(x_0)

                    let s_km1 = x_k - x_km1
                        y_km1 = grad_fx_k - grad_fx_km1
                        rho_km1 = 1 / (y_km1 `L.dot` s_km1) -- Scalar
                        v_km1 = L.ident (length gradEval) - (rho_km1 `L.scale` y_km1 `L.outer` s_km1) -- Scaling can happen before outer
                        h_k = (L.tr (v_km1) L.<> h_km1 L.<> v_km1) + (rho_km1 `L.scale` s_km1 `L.outer` s_km1)
                        gradPreconditioned = h_k L.#> grad_fx_k
                    in (L.toList gradPreconditioned,
                        bfgsParams { lastState = Just x_k, lastGrad = Just grad_fx_k, invH = Just h_k })

              _ -> error "invalid BFGS state"

      LBFGS -> -- Approximate the inverse of the Hessian times the gradient
               -- Only using the last `m` gradient/state difference vectors, not building the full h_k matrix (Nocedal p226)
            let grad_prev = lastGrad bfgsParams
                x_prev = lastState bfgsParams
                ss_val = s_list bfgsParams
                ys_val = y_list bfgsParams
                km1 = numUnconstrSteps bfgsParams -- Our current step is k; the last step is km1 (k_minus_1)
                m = memSize bfgsParams
            in case (grad_prev, x_prev, ss_val, ys_val, km1) of

               -- Perform normal gradient descent on first step
               (Nothing, Nothing, [], [], 0) ->
                   -- Store x_k, grad f(x_k) so we can compute s_k, y_k on next step
                   let bfgsParams' = bfgsParams { lastState = Just x_k, lastGrad = Just grad_fx_k,
                                                  s_list = [], y_list = [], numUnconstrSteps = km1 + 1 } in
                   (gradEval, bfgsParams')

               (Just grad_fx_km1, Just x_km1, ss, ys, km1) ->
                   -- Compute s_{k-1} = x_k - x_{k-1} and y_{k-1} = (analogous with grads)
                   -- Unlike Nocedal, compute the difference vectors first instead of last (same result, just a loop rewrite)
                   -- Use the updated {s_i} and {y_i}. (If k < m, this reduces to normal BFGS, i.e. we use all the vectors so far)
                   let (s_km1, y_km1) = (x_k - x_km1, grad_fx_k - grad_fx_km1) -- Newest vectors added to front
                       (ss', ys') = (take m $ s_km1 : ss, take m $ y_km1 : ys) -- The limited-memory part: drop stale vectors
                       gradPreconditioned = lbfgs grad_fx_k ss' ys'
                       descentDirCheck = -gradPreconditioned `L.dot` grad_fx_k

                   -- Reset L-BFGS if the result is not a descent direction, and use steepest descent direction
                   -- https://github.com/JuliaNLSolvers/Optim.jl/issues/143 https://github.com/JuliaNLSolvers/Optim.jl/pull/144

                   in if descentDirCheck < epsd
                      then (L.toList gradPreconditioned, bfgsParams { lastState = Just x_k, lastGrad = Just grad_fx_k,
                                                                      s_list = ss', y_list = ys', numUnconstrSteps = km1 + 1 })
                      else tro "L-BFGS did not find a descent direction. Resetting correction vectors." $
                           (gradEval, bfgsParams { lastState = Just x_k, lastGrad = Just grad_fx_k,
                                                   s_list = [], y_list = [], numUnconstrSteps = km1 + 1 })

                   -- TODO: check the curvature condition y_k^T s_k > 0 (8.7) (Nocedal 201)
                   -- https://github.com/JuliaNLSolvers/Optim.jl/issues/26

               _ -> error "invalid L-BFGS state"

type LVector = L.Vector L.R
type LMatrix = L.Matrix L.R

-- See Nocedal p225
-- expects ys and ss to be in order from most recent to oldest (k-1 ... k-m)
-- expects length ys == length ss, length ys > 0, length ys > 0
lbfgs :: LVector -> [LVector] -> [LVector] -> LVector
lbfgs grad_fx_k ss ys =
    let rhos = map calculate_rho $ zip ss ys -- The length of any list should be the number of stored vectors
        q_k = grad_fx_k
        (q_k_minus_m, alphas) = foldl' pull_q_back (q_k, []) (zip3 rhos ss ys) -- backward: for i = k-1 ... k-m
        -- Note the order of alphas will be from k-m through k-1 for the push_r_forward loop
        h_0_k = estimate_hess (head ys) (head ss)
        r_k_minus_m = h_0_k L.#> q_k_minus_m
        r_k = foldl' push_r_forward r_k_minus_m (zip (zip (reverse rhos) alphas) (zip (reverse ss) (reverse ys)))
        -- forward: for i = k-m .. k-1 (TODO: optimize out the reverses)
    in r_k -- is H_k * grad f(x_k)

           where calculate_rho :: (LVector, LVector) -> L.R
                 calculate_rho (s, y) = 1 / ((y `L.dot` s) + epsd)

                 pull_q_back :: (LVector, [L.R]) -> (L.R, LVector, LVector) -> (LVector, [L.R]) -- from i+1 to i
                 pull_q_back (q_i_plus_1, alphas) (rho_i, s_i, y_i) =
                        let alpha_i = rho_i * (s_i `L.dot` q_i_plus_1) -- scalar
                            q_i = q_i_plus_1 - (alpha_i `L.scale` y_i) -- scalar * vector
                        in (q_i, alpha_i : alphas) -- Note the order of alphas

                 -- Scale I by an estimate of the size of the Hessian along the most recent search direction (Nocedal p226)
                 estimate_hess :: LVector -> LVector -> LMatrix
                 estimate_hess y_km1 s_km1 =
                        let gamma_k = (s_km1 `L.dot` y_km1) / ((y_km1 `L.dot` y_km1) + epsd)
                        in gamma_k `L.scale` (L.ident (L.size y_km1))

                 push_r_forward :: LVector -> ((L.R, L.R), (LVector, LVector)) -> LVector -- from i to i+1
                 push_r_forward r_i ((rho_i, alpha_i), (s_i, y_i)) =
                        let beta_i = rho_i * (y_i `L.dot` r_i) -- scalar
                            r_i_plus_1 = r_i + (alpha_i - beta_i) `L.scale` s_i -- scalar * vector
                        in r_i_plus_1

estimateGradient :: ObjFn1 a -> [Float] -> [Float]
estimateGradient f state =
      let len = length state
          h = 0.001 -- Choice of h really matters!!! This is not an accurate estimate.
          fx = f state
      -- time is O(|state|^2)
          dfx i = ((f $ replace i ((state !! i) + h) state) - fx) / h
          dfxs = map dfx [0..(len-1)]
      in dfxs
      where replace pos newVal list = take pos list ++ newVal : drop (pos+1) list

-- Given the objective function, gradient function, timestep, and current state,
-- return the timestep (found via line search) and evaluated gradient at the current state.
-- the autodiff library requires that objective functions be polymorphic with Floating a
timeAndGrad :: OptConfig -> Params -> ObjFn1 a -> [Float] -> (Float, [Float], [Float], BfgsParams)
timeAndGrad config params f state =
            let gradEval = if useAutodiff then gradF state else estimateGradient f state
                (gradToUse_d, bfgs') = gradP config (bfgsInfo params) (map r2f gradEval :: [Double]) f state
                gradToUse = map r2f gradToUse_d
                -- Use line search to find a good timestep. If we use Newton's method, the descent direction uses the preconditioned gradient.
                descentDir = negL $ gradToUse
                timestep = let resT = if useLineSearch && useAutodiff
                                      then awLineSearch config f (duf descentDir) descentDir state
                                      else constT in -- hardcoded timestep
                           if isNaN resT then error "returned timestep is NaN" else resT

            in tr "timeAndGrad: " (timestep, gradEval, gradToUse, bfgs')

            where gradF :: GradFn a
                  gradF = appGrad f

                  -- directional derivative of f at x in the direction of u (descent direction, which may not have unit norm)
                  duf :: (Autofloat a) => [a] -> [a] -> a
                  duf u x = u `dotL` gradF x

linesearch_max :: Int
linesearch_max = 100 -- TODO what's a reasonable limit (if any)?

-- Implements Armijo-Wolfe line search as specified in Keenan's notes, converges on nonconvex fns as well
-- based off Lewis & Overton, "Nonsmooth optimization via quasi-Newton methods
-- duf = D_u(f), the directional derivative of f at descent direction u
-- D_u(x) = <gradF(x), u>. If u = -gradF(x) (as it is here), then D_u(x) = -||gradF(x)||^2

awLineSearch :: OptConfig -> ObjFn1 a -> ObjFn1 a -> [Float] -> [Float] -> Float
awLineSearch config f duf descentDir x0 =

     let (a0, b0, t0) = (0, infinity, 1) -- Unit step length should be tried first for quasi-Newton methods. (Nocedal 201)

         -- drop while a&w are not satisfied OR the interval is large enough
         (numUpdates, (af, bf, tf)) = head $ dropWhile intervalOK_or_notArmijoAndWolfe
                                      $ zip [0..] $ iterate update (a0, b0, t0)

     in tro ("Line search # updates: " ++ show numUpdates) $ tf

          where update :: (Float, Float, Float) -> (Float, Float, Float)
                update (a, b, t) =
                       let (a', b', sat) | not $ armijo t    = trl "not armijo" (a, t, False)
                                         | not $ wolfe t     = trl "not wolfe" (t, b, False)
                                         | otherwise         = (a, b, True) in
                       if sat then (a, b, t) -- if armijo and wolfe, then we use (a, b, t) as-is
                       else if b' < infinity then trl "b' < infinity" (a', b', (a' + b') / 2)
                       else trl "b' = infinity" (a', b', 2 * a')

                intervalOK_or_notArmijoAndWolfe :: (Int, (Float, Float, Float)) -> Bool
                intervalOK_or_notArmijoAndWolfe (numUpdates, (a, b, t)) =
                      not $
                      if armijo t && wolfe t then
                           trl ("stop: both sat. |descentDir at x0| = " ++ show (norm descentDir)) True
                      else if abs (b - a) < minInterval then
                           trl ("stop: interval too small. |descentDir at x0| = " ++ show (norm descentDir)) True
                      else if numUpdates > linesearch_max then
                           trl ("stop: number of line search updates exceeded max") True
                      else False

                armijo :: Float -> Bool
                armijo t = (f (x0 +. t *. descentDir)) <= (fAtx0 + c1 * t * dufAtx0)

                weakWolfe :: Float -> Bool -- Better for nonsmooth functions
                weakWolfe t = (duf (x0 +. t *. descentDir)) >= (c2 * dufAtx0)

                strongWolfe :: Float -> Bool
                strongWolfe t = (abs (duf (x0 +. t *. descentDir))) <= (c2 * abs dufAtx0)

                -- Descent direction (at x0) must have a negative dot product with the gradient of f at x0.
                dufAtx0 :: Float -- TODO: this is redundant, cache it
                dufAtx0 = let res = duf x0 in
                        trl ("<grad f(x0), descent direction at x0>: " ++ show res) $
                        if res > 0 then tro "WARNING: descent direction doesn't satisfy condition" res else res

                fAtx0 :: Float
                fAtx0 = f x0

                minInterval :: Float -- stop if the interval gets too small; might not terminate
                minInterval = if intervalMin then 10 ** (-10) else 0
                -- TODO: the line search is very sensitive to this parameter. Blows up with 10**(-5). Why?

                wolfe :: Float -> Bool
                wolfe = weakWolfe

max_stars_repo_name: daodaoliang/penrose

max_issues_repo_name: daodaoliang/penrose

max_forks_repo_name: daodaoliang/penrose


text:
module Utils where

import Data.Function (on)
import Data.List (sortBy)
import System.Random.MWC (GenIO,uniform)
import Control.Monad.Primitive (PrimMonad, PrimState)
import Numeric.LinearAlgebra.HMatrix (Vector)
import qualified Data.Vector.Storable as V

multinomial :: Vector Double  -- ^ Probabilities of each of the p different outcomes. These should sum to 1.
            -> GenIO
            -> IO (Vector Int)
multinomial pvals gen
    | V.sum pvals > 1.01 =
        error "Utils.multinomial: sum of probabilities should sum to 1"
    | otherwise = do
        y <- uniform gen
        let (_,sample,_) = prepvals y
        return sample
    where prepvals :: Double -> (Double, Vector Int, Bool)
          prepvals x = foldr
                  (\(j,p) (s,v,z) ->
                        if z
                            then (s, v, z)
                            else if s + p >= x
                                     then ( s
                                          , v V.//
                                            [(j, 1)]
                                          , True)
                                     else (s + p, v, z))
                  (0.0, V.replicate k 0, False)
                  (sortBy
                       (compare `on` snd)
                       (zip
                            [0 ..]
                            (V.toList pvals)))
          k = V.length pvals

ext: hs

max_stars_repo_name: DbIHbKA/mnistvb

max_issues_repo_name: DbIHbKA/mnistvb

max_forks_repo_name: DbIHbKA/mnistvb


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

module Feldspar.Core.Frontend.Floating where

import qualified Prelude
import Prelude (Float,Double)
import Data.Complex

import Feldspar.Core.Constructs
import Feldspar.Core.Constructs.Floating
import Feldspar.Core.Frontend.Fractional

-- Make new class, with "Data" in all the types

infixr 8 **

class (Fraction a, Prelude.Floating a) => Floating a where
  pi        :: Data a
  pi        =  sugarSymF Pi
  exp       :: Data a -> Data a
  exp       =  sugarSymF Exp
  sqrt      :: Data a -> Data a
  sqrt      =  sugarSymF Sqrt
  log       :: Data a -> Data a
  log       =  sugarSymF Log
  (**)      :: Data a -> Data a -> Data a
  (**)      =  sugarSymF Pow
  logBase   :: Data a -> Data a -> Data a
  logBase   =  sugarSymF LogBase
  sin       :: Data a -> Data a
  sin       =  sugarSymF Sin
  tan       :: Data a -> Data a
  tan       =  sugarSymF Tan
  cos       :: Data a -> Data a
  cos       =  sugarSymF Cos
  asin      :: Data a -> Data a
  asin      =  sugarSymF Asin
  atan      :: Data a -> Data a
  atan      =  sugarSymF Atan
  acos      :: Data a -> Data a
  acos      =  sugarSymF Acos
  sinh      :: Data a -> Data a
  sinh      =  sugarSymF Sinh
  tanh      :: Data a -> Data a
  tanh      =  sugarSymF Tanh
  cosh      :: Data a -> Data a
  cosh      =  sugarSymF Cosh
  asinh     :: Data a -> Data a
  asinh     =  sugarSymF Asinh
  atanh     :: Data a -> Data a
  atanh     =  sugarSymF Atanh
  acosh     :: Data a -> Data a
  acosh     =  sugarSymF Acosh

instance Floating Float
instance Floating Double

instance (Fraction a, Prelude.RealFloat a) => Floating (Complex a)

max_stars_repo_name: kffaxen/feldspar-language

max_issues_repo_name: kffaxen/feldspar-language

max_forks_repo_name: kffaxen/feldspar-language


text:
{-# LANGUAGE ScopedTypeVariables #-}

module FM (FMOptions(..), run, parseOpts) where

import           Control.Monad
import           Control.Monad.State.Lazy (StateT(runStateT), MonadState (get, put))
import           Data.Complex
import qualified Data.Vector.Storable as V
import qualified Filter as F
import           IQ
import           Options.Applicative
import           Pipes
import           Render
import qualified Sound.File.Sndfile as SF
import qualified Sound.File.Sndfile.Buffer.Vector as SV
import           System.IO
import           ZeroStuff

data FMOptions = FMOptions
  { inputFile :: String
  , outputFile :: String
  , deviationRatio :: Double
  , upsample :: Int
  , filterOpts :: F.FilterOptions }

parseOpts :: Parser FMOptions
parseOpts = FMOptions
         <$> strOption
             ( long "in-file"
            <> short 'i'
            <> metavar "INPUT_FILE"
            <> help "Input sound file" )
         <*> strOption
             ( long "out-file"
            <> short 'o'
            <> metavar "OUTPUT_FILE"
            <> help "Output IQ file")
         <*> option auto
             ( long "deviation-ratio"
            <> short 'd'
            <> value 1
            <> showDefault
            <> metavar "DEV_RATIO"
            <> help "The raio of max_devation over max_mod_frequency.")
         <*> option auto
             ( long "upsample-factor"
            <> short 'u'
            <> value 10
            <> showDefault
            <> metavar "UPSAMPLE_FAC"
            <> help "Upsample factor, outSR = inSR * UPSAMPLE_FAC")
         <*> F.parseOpts

-- clamp phase to [-pi, pi]
clampPhase :: Double -> Double
clampPhase phase
  | phase > pi = clampPhase $ phase - (2 * pi)
  | phase < -pi = clampPhase $ phase + (2 * pi)
  | otherwise =  phase

modulateFM :: MonadIO  m => Double -> Pipe Float IQ m ()
modulateFM devRatio = void $ flip runStateT (0 :: Double) $ forever $ do
  sample <- lift await
  phase <- get
  let phase' = clampPhase $ phase + realToFrac sample * devRatio
  lift . yield $ realToFrac (cos phase') :+ realToFrac (sin phase')
  put phase'

yieldEvery :: Monad m => Int -> Pipe a a m ()
yieldEvery n = do
  replicateM_ (n - 1) await
  await >>= yield
  yieldEvery n

run :: FMOptions -> IO ()
run opts = do
  (info, Just (aSamps :: SV.Buffer Float)) <- SF.readFile $ inputFile opts

  let sampsVec = SV.fromBuffer aSamps
      upsampleRate = upsample opts * SF.samplerate info
      filterKernel = F.lowPass upsampleRate $ filterOpts opts

  putStr "Processing samples... "
  hFlush stdout

  withBinaryFile (outputFile opts) WriteMode  $ \f ->
    runEffect $ V.mapM_ yield sampsVec
            >-> yieldEvery (SF.channels info)
            >-> modulateFM (deviationRatio opts)
            >-> zeroStuff (upsample opts)
            >-> F.convolve filterKernel
            >-> cfileSink f

  putStrLn "Done."

  putStrLn $  "CFile sample rate: " ++ show  upsampleRate ++ "Hz"

max_stars_repo_name: hexagonal-sun/ayeQ

max_issues_repo_name: hexagonal-sun/ayeQ

max_forks_repo_name: hexagonal-sun/ayeQ


text:
import           Numeric.LinearAlgebra
import           System.Random ( mkStdGen )
import           Learning ( nrmse )

import           RC.NTC as RC

-- Get training data
takePast :: Element a => Int -> Matrix a -> Matrix a
takePast horizon xs = xs ?? (All, Take (len - horizon))
  where
    len = cols xs

-- Get target data
takeFuture :: Element a => Int -> Matrix a -> Matrix a
takeFuture horizon = (?? (All, Drop horizon))

-- Manipulate the data matrix: horizontal split
splitAtRatio :: Element a => Double -> Matrix a -> (Matrix a, Matrix a)
splitAtRatio ratio m = splitAt' spl m
  where
    total = cols m
    spl = round $ ratio * fromIntegral total

splitAt' :: Element a => Int -> Matrix a -> (Matrix a, Matrix a)
splitAt' i m = (m ?? (All, Take i), m ?? (All, Drop i))

main :: IO ()
main = do
  -- Load and transpose time series to predict
  dta <- tr <$> loadMatrix "examples/data/mg.txt"

  -- Train on 50% of data
  let (train, validate) = splitAtRatio 0.50 dta

  -- Configure a new NTC network
  let p = RC.par0 { RC._inputWeightsRange = (0.1, 0.3) }
      g = mkStdGen 1111
      ntc = RC.new g p (1, 1000, 1)

  -- 20 steps ahead prediction horizon
  let horizon = 20

  let train' = takePast horizon train  -- Past series
      trainTeacher = takeFuture horizon train  -- Predicted series

  let forgetPts = 300  -- Washout

  -- Train
  case RC.learn ntc forgetPts train' trainTeacher of
    Left s -> error s
    Right ntc' -> do
      let target = (takeFuture horizon validate) ?? (All, Drop forgetPts)

      -- Predict
      case RC.predict ntc' forgetPts (takePast horizon validate) of
        Left s -> error s
        Right prediction -> do
          let tgt' = flatten target
              predic' = flatten prediction
              err = nrmse tgt' predic'

          putStrLn $ "Error: " ++ show err

          let result = (tr target) ||| (tr prediction)

          saveMatrix "examples/NTC/result.txt" "%g" result

max_stars_repo_name: masterdezign/rc

max_issues_repo_name: masterdezign/rc

max_forks_repo_name: masterdezign/rc


text:
{-# LANGUAGE BangPatterns #-}

module SimParse where

import Data.List
import NEC
import System.Process
import System.Random
import Data.Complex

vswr :: Double -> (Double, Double) -> Double
vswr target = vswr' (target, 0)

vswr' :: (Double, Double) -> (Double, Double) -> Double
vswr' (reTar, imTar) (reLoad, imLoad) = swr where
  z0 = reTar :+ imTar
  zL = reLoad :+ imLoad
  refRatio = (zL - z0) / (zL + z0)
  (ro :+ _) = abs $ refRatio
  swr = (1 + ro)  / (1 - ro)

bandwidthInterval :: Double -> Double -> [[String]] -> (Double, Double)
bandwidthInterval targetImp swrThre lnss = (lo, hi) where
  swrs = map (\l -> (frequency l , vswr targetImp (impedance l))) lnss
  below = takeWhile (\s -> (snd s) <= swrThre) $ dropWhile (\s -> (snd s) > swrThre) swrs
  lo = fst $ head below
  hi = fst $ last below

bandwidth t s = (\(lo, hi) -> hi - lo) . bandwidthInterval t s

gain :: [String] -> Double
gain lns = read gn where
  gn = head $ drop 4 $ words gainline
  gainline =  head $ drop 5 $ dropWhile (not . isSubsequenceOf "RADIATION PATTERN") lns

impedance :: [String] -> (Double, Double)
impedance lns = (zr, zi) where
  zr = read $ head line
  zi = read $ head $ tail line
  line = drop 6 $ words paramsLine
  paramsLine = head $ drop 3 $ dropWhile (not . isSubsequenceOf "ANTENNA INPUT PARAMETERS") lns

frequency :: [String] -> Double
frequency lns = read freq where
  freq = head $ drop 2 $ words $ head $ dropWhile (not . isSubsequenceOf "FREQUENCY") $ lns

fitness targetImpedance lns = a * gn - b * (abs $ targetImpedance - reZ) - c * (abs imZ) where
  a = 40
  b = 2
  c = 2
  gn = gain lns
  (reZ, imZ) = impedance lns

fitnessImepd targetImpedance lns = - b * (abs $ targetImpedance - reZ) - c * (abs imZ) where
  b = 2
  c = 2
  gn = gain lns
  (reZ, imZ) = impedance lns

fitness' targetImpedance lns = a * gn - b * ((targetImpedance - reZ)*(targetImpedance - reZ)) - c * (imZ * imZ) where
  a = 40
  b = 2
  c = 2
  gn = gain lns
  (reZ, imZ) = impedance lns

fitnessImepd' targetImpedance lns =  - b * ((targetImpedance - reZ)*(targetImpedance - reZ)) - c * (imZ * imZ) where
  b = 2
  c = 2
  gn = gain lns
  (reZ, imZ) = impedance lns

runSim sim sfreq efreq steps = go where
  splitRuns [] = []
  splitRuns xs = first : splitRuns rest where
    first = takeWhile (not . isSubsequenceOf "- FREQUENCY -") $ rest
    rest = tail $ (dropWhile (not . isSubsequenceOf "- FREQUENCY -")) xs

  go = do
    rn <- randomIO :: IO Int
    let fl = "/tmp/evalnec"++(show rn)
    writeFile (fl++".nec") $ printSim sfreq efreq steps False sim
    system $ "nec2c -i "++fl++".nec"
    !lns <- fmap lines $ readFile (fl++".out")
    system $ "rm "++fl++"*"
    return $ init $ splitRuns $ lns ++ ["- FREQUENCY -"]

evalSim sim freq fitnessF = go where
  --freq = 447.625
  go = do
    rn <- randomIO :: IO Int
    let fl = "/tmp/evalnec"++(show rn)
    writeFile (fl++".nec") $ printSim freq freq 1 False sim
    system $ "nec2c -i "++fl++".nec"
    !lns <- fmap lines $ readFile (fl++".out")
    system $ "rm "++fl++"*"
    --print $ gain lns
    --print $ impedance lns
    return $ fitnessF lns

max_stars_repo_name: bssstudio/NECGenHS

max_issues_repo_name: bssstudio/NECGenHS

max_forks_repo_name: bssstudio/NECGenHS


text:
module DataTypes where

import Control.Monad.Except
import Data.Array
import Data.Complex
import Data.IORef
import Data.Ratio
import GHC.IO.Handle
import ParserCombinators (ParseError)

type Env = IORef EnvFrame

type Bindings = [(String, IORef LispVal)]

data EnvFrame
  = Global { bindings :: Bindings }
  | Frame { parent :: Env
          , bindings :: Bindings }

type IOThrowsError = ExceptT LispError IO

type PrimitiveFunc = ([LispVal] -> ThrowsError LispVal)

type IOPrimitiveFunc = ([LispVal] -> IOThrowsError LispVal)

type ThrowsError = Either LispError

data LispError
  = NumArgs Arity
            Int
            [LispVal]
  | TypeMismatch String
                 LispVal
  | Parser ParseError
  | BadSpecialForm String
                   LispVal
  | NotFunction String
                String
  | UnboundVar String
               String
  | Default String

instance Show LispError where
  show = showError

data LispVal
  = Vector (Array Integer LispVal)
  | Atom String
  | List [LispVal]
  | DottedList [LispVal]
               LispVal
  | Integer Integer
  | Rational Rational
  | Float Double
  | Complex (Complex Double)
  | String String
  | Character Char
  | Bool Bool
  | IOFunc String
           IOPrimitiveFunc
  | Port Handle
  | PrimitiveFunc String
                  PrimitiveFunc
  | Func { name :: String
         , params :: [String]
         , vararg :: (Maybe String)
         , body :: [LispVal]
         , closure :: Env }
  | Void

instance Show LispVal where
  show = showVal

data Arity
  = Min Int
  | MinMax Int
           Int

unwordsList :: [LispVal] -> String
unwordsList = unwords . map showVal

showVal :: LispVal -> String
showVal (String contents) = "\"" ++ contents ++ "\""
showVal (Atom name') = name'
showVal (Integer contents) = show contents
showVal (Rational contents) =
  show (numerator contents) ++ "/" ++ show (denominator contents)
showVal (Float contents) = show contents
showVal (Complex contents) =
  show (realPart contents) ++
  (if imagPart contents < 0
     then ""
     else "+") ++
  show (imagPart contents) ++ "i"
showVal (Character contents) = show contents
showVal (Bool True) = "#t"
showVal (Bool False) = "#f"
showVal (List contents) = "(" ++ unwordsList contents ++ ")"
showVal (Vector contents) = "#(" ++ unwordsList (elems contents) ++ ")"
showVal (DottedList h t) = "(" ++ unwordsList h ++ " . " ++ showVal t ++ ")"
showVal (PrimitiveFunc name' _) = "#<procedure:" ++ name' ++ ">"
showVal Func {name = name'} = "#<procedure:" ++ name' ++ ">"
showVal (Port _) = "<IO port>"
showVal (IOFunc name' _) = "#<IO primitive:" ++ name' ++ ">"
showVal Void = ""

showError :: LispError -> String
showError (Default message) = message
showError (UnboundVar message varname) = message ++ ": " ++ varname
showError (BadSpecialForm message form) = message ++ ": " ++ show form
showError (NotFunction message func) = message ++ ": " ++ show func
showError (NumArgs expected found args) =
  let error' =
        case expected of
          Min min' ->
            "arity mismatch;\nthe expected number of arguments does not match the given number" ++
            "\nexpected: at least " ++ show min' ++ "\ngiven: " ++ show found
          MinMax min' max' ->
            "arity mismatch;\nthe expected number of arguments does not match the given number" ++
            "\nexpected: " ++
            (if min' == max'
               then show min'
               else "between " ++ show min' ++ " and " ++ show max') ++
            "\ngiven: " ++ show found
      argsError =
        case args of
          [] -> ""
          a -> "\narguments:\n" ++ unwordsList a
  in error' ++ argsError
showError (TypeMismatch expected found) =
  "Invalid type: expected " ++ expected ++ ", found " ++ show found
showError (Parser parseErr) = "Parse error at " ++ show parseErr

trapError :: (Show a, MonadError a m) => m String -> m String
trapError action = catchError action (return . show)

extractValue :: ThrowsError a -> a
extractValue (Right val) = val

max_stars_repo_name: ebenpack/hascheme

max_issues_repo_name: ebenpack/hascheme

max_forks_repo_name: ebenpack/hascheme


text:
{-# LANGUAGE FlexibleContexts #-}
{-|
Module      : $Header$
CopyRight   : (c) 8c6794b6, 2011-2013
License     : BSD3
Maintainer  : 8c6794b6@gmail.com
Stability   : experimental
Portability : portable

Synthesis methods

-}
module Spectrofy.Synth
  ( sinsyn
  , fftsyn
  ) where

import Data.Complex (Complex(..), realPart)
import Data.Word (Word8)
import Data.List (foldl1')

import Data.Array.Repa
    ((:.)(..), Array, All(..), D, DIM2, DIM3, Source, U, Z(..))
import qualified Data.Array.Repa as R
import qualified Data.Array.Repa.FFTW as R

-- --------------------------------------------------------------------------
--
-- Manual sinusoids
--

-- | Sums up sinusoids manually.
--
-- Frequency is taken from y axis, time from x axis.
-- Amplitude for each frequency is result of luminated value in each pixel.
--
sinsyn ::
    Source s (Word8, Word8, Word8)
    => Int -> Int -> Array s DIM2 (Word8, Word8, Word8) -> Array D DIM2 Double
sinsyn rate delta = squash . sumSins . toSins rate . to3D delta . rgbamp
{-# INLINE sinsyn #-}

rgbamp ::
    Source s (Word8, Word8, Word8)
    => Array s DIM2 (Word8, Word8, Word8) -> Array D DIM2 Double
rgbamp arr = {-# SCC "rgbamp" #-} R.traverse arr f g where
  f = id
  g h ix = 10 ** ((red*0.21 + green*0.71 + blue*0.07) / (255*3)) - 1 where
    (red', green', blue') = h ix
    red = fromIntegral red'
    green = fromIntegral green'
    blue = fromIntegral blue'
{-# INLINE rgbamp #-}

to3D :: Source s Double => Int -> Array s DIM2 Double -> Array D DIM3 Double
to3D n = R.extend (Z:.All:.All:.n)
{-# INLINE to3D #-}

toSins :: Source s Double => Int -> Array s DIM3 Double -> Array D DIM3 Double
toSins rate arr = {-# SCC "toSins" #-} R.traverse arr id f where
  _:.y:._:._ = R.extent arr
  f g ix@(_:.i:._:.k)
    | amp <= 0  = 0
    | otherwise = amp * sin (frq * k' * 2 * pi / rate') / y'
    where
      amp = g ix
      frq = (rate'/2) - (((y'-i'+1)/y') * (rate'/2))
      rate' = fromIntegral rate
      i' = fromIntegral i
      y' = fromIntegral y
      k' = fromIntegral k
{-# INLINE toSins #-}

sumSins :: (Source s Double) => Array s DIM3 Double -> Array U DIM2 Double
sumSins arr = {-# SCC "sumSins" #-} R.sumS (R.backpermute sh' f arr) where
  f (_:.i:.j:.k) = Z:.k:.j:.i
  sh' = Z:.z:.y:.x
  (_:.x:.y:.z) = R.extent arr
{-# INLINE sumSins #-}

squash :: Source s Double => Array s DIM2 Double -> Array D DIM2 Double
squash arr = {-# SCC squash #-} R.backpermute sh' f arr where
  f (_:._:.j) = Z:.j `mod` x:.j `div` x
  (_:.x:.y) = R.extent arr
  sh' = Z:.1:.(x*y)
{-# INLINE squash #-}

-- --------------------------------------------------------------------------
--
-- Inversed FFT
--

-- | Performs 1-dimensional inverse fft on each column, then concatenate results.
--
-- RGB values are used for initial phases and magnitude.
--
fftsyn :: Int -> Array D DIM2 (Word8, Word8, Word8) -> Array D DIM2 Double
fftsyn fsz img =
  {-# SCC "fftsyn" #-}
  foldl1' R.append $ map (fftslice fsz img') [0..n-1] where
    img' = rgbamp img
    _:._:.n = R.extent img
{-# INLINE fftsyn #-}

fftslice :: Int -> Array D DIM2 Double -> Int -> Array D DIM2 Double
fftslice wsize arr n =
  {-# SCC "fs_reshape" #-}
  R.reshape (Z :. 1 :. wsize :: DIM2) $
  {-# SCC "fs_map" #-}
  R.map realPart $
  {-# SCC "fs_ifft" #-}
  R.ifft $
  {-# SCC "fs_computeS" #-}
  R.computeS $
  {-# SCC "fs_slice" #-}
  R.slice (R.map (\x -> x :+ x) $ grow2d wsize (zpad arr)) (Z:.All:.n)
{-# INLINE fftslice #-}

grow2d :: Int -> Array D DIM2 Double -> Array D DIM2 Double
grow2d n arr = {-# SCC "grow2d" #-} R.traverse arr f g where
  f _ = Z :. n :. y
  _:.x:.y = R.extent arr
  g h (_:.i:.j)
    | i' >= 4   = k/4
    | otherwise = h (Z:.i':.j)
    where
      k  = h (Z:.i'-3:.j) + h (Z:.i'-2:.j) + h (Z:.i'-1:.j) + h (Z:.i':.j)
      i' = fst $ properFraction $ (x' * fromIntegral i / n')
      n' = fromIntegral n :: Double
      x' = fromIntegral x
{-# INLINE grow2d #-}

zpad :: Array D DIM2 Double -> Array D DIM2 Double
zpad arr = {-# SCC "zpad" #-} R.reshape sh' arr' where
  _:.x:.y = R.extent arr
  len = x * y
  zeros = R.fromFunction zsh (const 0)
  zsh = Z :. len
  sh' = Z :. (2*x) :. y
  arr' = R.append (R.reshape zsh arr) zeros
{-# INLINE zpad #-}


text:
module IFS.Base
  ( IFS(..)
  , Transformation(..)
  , iter
  , chaosIter
  , sierpinskiTIFS
  , crossFrac
  , cantorSetIFS
  , squareLeafIFS
  , infTriangles
  , infTriangles2
  , infTriangles3
  , kochSnowFlakeIFS
  , differentFlakeIFS
  , differentFlake2IFS
  )
  where

import Numeric.LinearAlgebra.Data
import Numeric.LinearAlgebra.HMatrix
import System.Random
import Control.Applicative

type Transformation a = Vector a -> Vector a

type IFS = [Transformation R]

iter :: IFS -> Vector R -> [Vector R]
iter ifs vec = ifs <*> pure vec

chaosIter :: IFS -> StdGen -> Vector R -> (Vector R, StdGen)
chaosIter ifs g v = (randSim v, g')
  where
    (index, g') = randomR (0, length ifs - 1) g
    randSim = ifs !! index

rotate2D :: Double -> Matrix R
rotate2D theta = (2><2) [cos theta', -sin theta',
                         sin theta', cos theta']
  where theta' = theta * pi / 180

sierpinskiTIFS :: IFS
sierpinskiTIFS =
  [ \v -> 1/2 * v
  , \v -> 1/2 * v + vector [1/2, 0]
  , \v -> 1/2 * v + vector [1/4, sqrt 3 / 4]
  ]

cantorSetIFS :: IFS
cantorSetIFS =
  [ \v -> 1/3 * v
  , \v -> 1/3 * v + vector [2/3, 0]
  ]

squareLeafIFS :: IFS
squareLeafIFS =
  [ \v -> 2/3 * v
  , \v -> 1/3 * v + vector [0, 2/3]
  , \v -> 1/3 * v + vector [2/3, 0]
  , \v -> 1/3 * v + vector [2/3, 2/3]
  ]

infTriangles :: IFS
infTriangles =
  [ \v -> 1/2 - 1/2 * v
  , \v -> 1/2 * v + vector [1/2, 0]
  , \v -> 1/2 * v + vector [0, 1/2]
  ]

infTriangles2 :: IFS
infTriangles2 =
  [ \v -> 1/2 + vector [-1/2, 1/2] * v
  , \v -> 1/2 * v + vector [1/2, 0]
  , \v -> 1/2 * v + vector [0, 1/2]
  ]

infTriangles3 :: IFS
infTriangles3 =
  [ \v -> vector [1/2, 0] + vector [-1/2, 1/2] * v
  , \v -> 1/2 * v + vector [1/2, 0]
  , \v -> 1/2 * v + vector [0, 1/2]
  ]

crossFrac :: IFS
crossFrac =
  [ \v -> 1/4 * v
  , \v -> 1/4 * v + vector [3/4, 0]
  , \v -> 1/4 * v + vector [0, 3/4]
  , \v -> 1/4 * v + vector [3/4, 3/4]
  , \v -> 1/2 * v + vector [1/4, 1/4]
  ]

kochSnowFlakeIFS :: IFS
kochSnowFlakeIFS =
  [ \v -> 1/3 * v
  , \v -> 1/3 * v + vector [2/3, 0]
  , \v -> rotate2D (60) `app` (1/3 * v) + vector [1/3, 0]
  , \v -> rotate2D (-60) `app` (1/3 * v) + vector [1/2, sqrt 3 / 6]
  ]

differentFlakeIFS :: IFS
differentFlakeIFS =
  [ \v -> 1/3 * v + vector [1 / 6, sqrt 3]
  , \v -> 1/3 * v + vector [1 / 6, sqrt 3 / 9]
  , \v -> 1/3 * v + vector [8 / 6, sqrt 3 / 9]
  , \v -> 1/3 * v + vector [8 / 6, sqrt 3]
  , \v -> 2/3 * v + vector [1 / 3, 1/3 * sqrt 3]
  ]

differentFlake2IFS :: IFS
differentFlake2IFS =
  [ \v -> 1/2 * v + vector [1/2, 0]
  , \v -> 1/4 * v + vector [3/4, - sqrt 3 / 4]
  , \v -> 1/4 * v + vector [3/4, sqrt 3 / 2]
  , \v -> 1/4 * v + vector [1/3 - sqrt 3 / 4, - sqrt 3 / 24]
  , \v -> 1/4 * v + vector [1/3 - sqrt 3 / 4, sqrt 3 / 3 - sqrt 3 / 24]
  , \v -> 1/4 * v + vector [13/12, - sqrt 3 / 24]
  , \v -> 1/4 * v + vector [13/12, sqrt 3 / 3 - sqrt 3 / 24]
  ]


text:
module Main where
import Data.Word
import Data.List.Split
import Data.Scientific as Scientific
import Data.Complex
import Control.Monad
import Codec.Picture
import Data.List
import Data.Time.Clock
import Data.Time.Calendar
import Graphics.X11.Xlib
import System.Exit (exitWith, ExitCode(..))
import Control.Concurrent
import Data.Bits
import GHC.Conc (numCapabilities)
import Graphical
import Mand
import GenerateSet
import GenerateInParallel

main :: IO ()
main = do
  taking_input <- newMVar True
 -- writePng ("benchmark.png") $ generateImage (\x y -> generatingFunction  0.0005 (-2) (2) (2) x y) 8000 8000
  generateInParallel
 -- forever $ loop taking_input

loop taking_input = do
  --print $ head $ reverse $ general_list_R mand_iterationA ((scientific 24 (-2)) :+ scientific 2 (-1))
  --print $ head $ reverse $ general_list_D mand_iteration (0.24 :+ (0.2))

  --writePng ("benchmark.png") $ generateImage (\x y -> f 0.0005 (-2) (2) (2) x y) 8000 8000
  input <- promptForImput taking_input "set for generating mandlebrot set, mov for the movement of the mandelbrot, grp for graphical, rang for range in a gif, sci for using scientific data type "
  case input of
    "set" -> generateSet taking_input
    "mov" -> generateMovement
    "grp" -> forkIO (startGraphical taking_input)>> return ()
    "rang" -> range
    "sci" -> generateSciSet
    "series" -> generateSeriesOfSets
    "test" -> tester
    "acc" -> generateAccurateSet
    "par" -> generateInParallel

generateAccurateSet = print "yet again not finishe"

tester = do
  writePng ("benchmark.png") $ generateImage (\x y -> generatingFunction  0.0005 (-2) (2) (2) x y) 8000 8000
--stepf srf sif realend x y
generateSeriesOfSets = print "in complete"

generateSciSet = do
  print "nothign yet"
-- print $ mi (makeSciComplex "1e-1" "2e-1") $ makeSciComplex "-9e-1" "4e-2"

--mi :: CS -> CS -> CS
--mi c z = c + z*z

makeSciComplex :: String -> String -> CS
makeSciComplex x y = (read x :: S) :+ (read x :: S)

range = do
  print "in complete"

generateMovement = do
  let list = [fromIntegral x/10:+ fromIntegral y/10 |x <-[(-400)..399], y <- [(-400)..399]]
  let zlist = foldl (\acc x -> (0:+0):acc) [] [0..639999]
  let list' = map' mand_iteration list list
  --print list
  --writePng "movv5_i1.png" $ generateImage (genMovImg $ map roundComplex2 list) 399 399
  --writePng "movv5_i2.png" $ generateImage (genMovImg $ map roundComplex2 list') 399 399
  --let list'1 = filterC $ map' mand_iteration list list'
  --writePng "movv5_i3.png" $ generateImage (genMovImg $ map roundComplex2 list'1) 399 399
  --let list'2 = filterC $ map' mand_iteration list list'1
  --writePng "movv5_i4.png" $ generateImage (genMovImg $ map roundComplex2 list'2) 399 399
  --let list'3 = filterC $ map' mand_iteration list list'2
  --writePng "movv5_i5.png" $ generateImage (genMovImg $ map roundComplex2 list'3) 399 399
  --let list'4 = filterC $ map' mand_iteration list list'3

  mapM_ (curry' writePng) $ reverse $ snd $ mapAccumR nextMov (list, zlist) [0..10]
  print "complete"

curry' :: (a -> b -> c) -> (a, b) ->  c
curry' f (a,b) = f a b

nextMov (last, original) x = ((nextl, nextOrig), ("movv10_i"++(show x)++".png",generateImage (genMovImg $ map roundComplex2 nextl) 3999 3999))
  where
    (nextl, nextOrig) = filter2 (map' mand_iteration original last ) original removeLargeC

removeLargeC :: (Num a, Ord a )=> Complex a -> Bool
removeLargeC (a :+ b) = if abs a >2 || abs b > 2 then False else True

filterC = filter removeLargeC

filter2 :: [a] -> [b] -> (a -> Bool) -> ([a],[b])
filter2 x y f = (fst ziped,snd ziped)
  where
    ziped = foldr (\(a,b) (ac,bc) -> if f a then (a:ac,b:bc) else (ac,bc)) ([],[]) $ zip x y

genMovImg arr x y = if elem (fromIntegral (x-2000) /1000 :+ fromIntegral (y-2000)/1000) arr then PixelRGB8 0 0 0 else PixelRGB8 255 255 255

map' :: (a -> a -> b) -> [a] -> [a] -> [b]
map' _ [] [] = []
map' f (x:xs) (y:ys) = (f x y):(map' f  xs ys)

g :: [Bool]
g = foldr (\x acc -> if mod x 10 == 0 then True:acc else False:acc) [] [0..399]

h :: [[Bool]]
h = foldr (\x acc -> if mod x 10 == 0 then g:acc else i:acc) [] [0..399]

i :: [Bool]
i = foldr (\_ acc -> False:acc) [] [0..399]

falseArr = foldr (\_ acc -> i:acc) [] [0..399]

toComplex :: RealFloat a => [[Bool]] -> [Complex a]
toComplex arr = foldr (\x acc1 -> (foldr (\y acc2 -> if arr !! x !! y then (fromIntegral (x-200)/100 :+ fromIntegral (y-200)/100):acc2 else acc2) [] [0..((length $ arr !! 0)-1)])++acc1) [] [0..((length arr) -1)]

toArr :: RealFloat a => [Complex a] -> [[Bool]]
toArr arr = f arr falseArr
  where
    rarr = map roundComplex2 arr
    f :: RealFloat a=> [Complex a] -> [[Bool]]-> [[Bool]]
    f [] fa  = fa
    f (x:xs) fa = f xs $ modifyarr fa  newx newy
      where
        newx = round $ 1000*(realPart x) +2000
        newy = round $ 1000*(imagPart x) +2000

modifyarr :: [[Bool]] -> Int -> Int -> [[Bool]]
modifyarr arr x y = foldr (\a acc -> (foldr (\b acc2 -> if x == a && b == y then True:acc2 else (arr !! x !! y):acc2) [] [0..((length $ arr !! 0)-1)]):acc) [] [0..((length arr)-1)]

roundComplex2 :: RealFloat a => Complex a -> Complex a
roundComplex2 (a :+ b) = (fromIntegral (round $ a*1000)/1000 :+ fromIntegral (round $ b*1000)/1000)

max_stars_repo_name: AlwinHughes/fractal

max_issues_repo_name: AlwinHughes/fractal

max_forks_repo_name: AlwinHughes/fractal


text:
{-# OPTIONS_GHC -fplugin GHC.TypeLits.KnownNat.Solver #-}
{-# OPTIONS_GHC -fno-warn-missing-signatures #-}
{-# OPTIONS_GHC -fno-warn-incomplete-patterns #-}

{-# LANGUAGE CPP                 #-}
{-# LANGUAGE ConstraintKinds     #-}
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE KindSignatures      #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}
{-# LANGUAGE TypeApplications    #-}
{-# LANGUAGE TypeOperators       #-}

module Test.Grenade.Layers.Mul where

import           System.Random.MWC

import           Data.Proxy

import           GHC.TypeLits

import           Grenade.Core
import           Grenade.Layers.Mul
import qualified Numeric.LinearAlgebra        as LA
import           Numeric.LinearAlgebra.Static (L, R)
import qualified Numeric.LinearAlgebra.Static as H

import           Hedgehog
import qualified Hedgehog.Range               as Range

import           Test.Hedgehog.Compat
import           Test.Hedgehog.Hmatrix

import           Data.Either
import           Data.Serialize

prop_mul_scalar_one_does_nothing = property $ do
  height   :: Int <- forAll $ choose 2 100
  width    :: Int <- forAll $ choose 2 100
  channels :: Int <- forAll $ choose 2 100

  case (someNatVal (fromIntegral height), someNatVal (fromIntegral width), someNatVal (fromIntegral channels)) of
    (Just (SomeNat (Proxy :: Proxy h)), Just (SomeNat (Proxy :: Proxy w)), Just (SomeNat (Proxy :: Proxy c))) -> do
      input :: S ('D3 h w c) <- forAll genOfShape

      let layer    = initMul :: Mul 1 1 1
          S3D out  = snd $ runForwards layer input :: S ('D3 h w c)
          S3D inp' = input :: S ('D3 h w c)

      H.extract inp' === H.extract out

prop_mul_random_scalar_as_expected = property $ do
  height   :: Int <- forAll $ choose 2 100
  width    :: Int <- forAll $ choose 2 100
  channels :: Int <- forAll $ choose 2 100

  scalar :: R 1 <- forAll randomVector
  let scale = H.extract scalar LA.! 0

  case (someNatVal (fromIntegral height), someNatVal (fromIntegral width), someNatVal (fromIntegral channels)) of
    (Just (SomeNat (Proxy :: Proxy h)), Just (SomeNat (Proxy :: Proxy w)), Just (SomeNat (Proxy :: Proxy c))) -> do
      input :: S ('D3 h w c) <- forAll genOfShape

      let layer    = Mul scalar :: Mul 1 1 1
          S3D out  = snd $ runForwards layer input :: S ('D3 h w c)
          inp' = (\(S3D x) -> H.dmmap (* scale) x) input :: L (h * c) w

      H.extract inp' === H.extract out

prop_mul_has_show = withTests 1 $ property $ do
  gen <- evalIO create
  mul :: Mul 1 1 1 <- evalIO $ createRandomWith UniformInit gen
  show mul `seq` success

prop_mul_can_update = withTests 1 $ property $ do
  gen <- evalIO create
  mul :: Mul 1 1 1 <- evalIO $ createRandomWith UniformInit gen
  runUpdate defSGD mul () `seq` success
  runUpdate defAdam mul () `seq` success

prop_mul_can_be_used_with_batch = withTests 1 $ property $ do
  reduceGradient @(Mul 1 1 1) [()] `seq` success

prop_mul_can_be_serialized = withTests 1 $ property $ do
  s <- forAll $ genRealNum $ Range.constant 1 10
  let mul :: Mul 1 1 1 = Mul (H.fromList [s])
      bs = encode mul
      dec = decode bs :: Either String (Mul 1 1 1)
  assert $ isRight dec
  let Right mul' = dec
      (Mul scalarMat)  = mul
      (Mul scalarMat') = mul'
  (H.extract scalarMat) === (H.extract scalarMat')

prop_mul_fails_serialize_with_incorrect_size = withTests 1 $ property $ do
  s <- forAll $ genRealNum $ Range.constant 1 10
  let mul :: Mul 1 1 1 = Mul (H.fromList [s])
  let bs = encode mul
  let dec = decode bs :: Either String (Mul 1 2 1)
  assert $ isLeft dec

tests :: IO Bool
tests = checkParallel $$(discover)

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
{-# LANGUAGE OverloadedStrings #-}

module Main where

import Numeric.LinearAlgebra

main :: IO ()
main = putStrLn $ show $ vector [1,2,3] * vector [3,0,-2]

max_stars_repo_name: tweag/economics

max_issues_repo_name: tweag/economics

max_forks_repo_name: tweag/economics


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TypeFamilies #-}

module Song.Grenade (
  SongSD
, song2SD
, song2Rn
, song2TupleRn
) where

import Song (Song(Song, _genre, _loudness, _tempo, _timeSignature, _key, _mode, _duration, _avgTimbre, _varTimbre))
import Grenade (S(S1D), Shape(D1))
import Numeric.LinearAlgebra.Static (R, konst, vector, (#), (&))
import Control.Arrow ((***), first)

type SongSD = (S ('D1 30), S ('D1 1))

song2SD :: (String -> Double) -> Song -> SongSD
song2SD tag = (S1D *** S1D) . song2Rn tag

song2Rn :: (String -> Double) -> Song -> (R 30, R 1)
song2Rn tag song = first (uncurry (#)) $ song2TupleRn tag song

song2TupleRn :: (String -> Double) -> Song -> ((R 3, R 27), R 1)
song2TupleRn tag Song{..} = ((discreteFeats, floatFeats), konst $ tag _genre)
  where
    discreteFeats = konst (fromIntegral _timeSignature)
                  & fromIntegral _key
                  & fromIntegral _mode

    floatFeats = (konst _loudness & _tempo & _duration)
               # (vector _avgTimbre::R 12)
               # (vector _varTimbre::R 12)

max_stars_repo_name: helq/haskell-binary-classification

max_issues_repo_name: helq/haskell-binary-classification

max_forks_repo_name: helq/haskell-binary-classification


text:
{-# LANGUAGE TemplateHaskell #-}

module School.Unit.Test.WeightDecay
( weightDecayTest ) where

import Conduit(Identity, liftIO)
import Numeric.LinearAlgebra ((><))
import School.TestUtils (diffCost, randomMatrix)
import School.Types.FloatEq ((~=))
import School.Types.Slinky (Slinky(..))
import School.Unit.CostFunction (CostFunction(..))
import School.Unit.UnitActivation (UnitActivation(..))
import School.Unit.UnitGradient (UnitGradient(..))
import School.Unit.WeightDecay
import School.Utils.LinearAlgebra (compareDoubleMatrix)
import Test.Tasty (TestTree)
import Test.Tasty.QuickCheck hiding ((><))
import Test.Tasty.TH
import Test.QuickCheck.Monadic (monadicIO, assert, pre)

eps :: Double
eps = 1e-5

prec :: Double
prec = 5e-2

weight :: Double -> CostFunction Double Identity
weight coeff = weightDecay coeff

prop_coeff_zero :: (Positive Int) -> (Positive Int) -> Property
prop_coeff_zero (Positive bSize) (Positive fSize) = monadicIO $ do
  inputMat <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation inputMat
  let cost = computeCost (weight 0) input SNil
  either (const (assert False))
         (\result -> assert $ result ~= 0)
         cost

prop_numerical_deriv :: (Positive Int) -> (Positive Int) -> Double -> Property
prop_numerical_deriv (Positive bSize) (Positive fSize) coeff = monadicIO $ do
  pre $ bSize < 23 && fSize < 23
  inputMat <- liftIO $ randomMatrix bSize fSize
  let input = BatchActivation inputMat
  let costFunc = weight coeff
  let deriv = derivCost costFunc input SNil
  let check = (bSize >< fSize) [ diffCost costFunc
                                          input
                                          eps
                                          (j, k)
                                          SNil
                               | j <- [0..bSize-1], k <- [0..fSize-1] ]
  either (const (assert False))
         (\(BatchGradient g) -> assert $ compareDoubleMatrix prec g check)
         deriv

weightDecayTest :: TestTree
weightDecayTest = $(testGroupGenerator)

max_stars_repo_name: jfulseca/School

max_issues_repo_name: jfulseca/School

max_forks_repo_name: jfulseca/School


text:
{-# LANGUAGE MultiParamTypeClasses #-}
-- | Defines intersection tests for a varity of geometry.
module Game.Physics.Collision ( Point
                              , point2D
                              , AABB
                              , aabb2D
                              , OBB
                              , obb2DCorners
                              , Intersectable(..)
                              ) where

import Prelewd hiding (toList, (<>))

import Impure

import Data.Tuple
import Storage.List
import Numeric.LinearAlgebra

-- | "Axis-aligned bounding box."
--
--   Use the smart constructor 'aabb2D' to build AABBs.
--
--   * Implementation details:
--
--   An AABB is an Nx2 matrix. X[0] is the lower limit on the x-th dimension
--   and X[1] is the upper limit on the x-th dimension.
--
--   The number of rows is the number of dimensions of the box. The number of
--   columns is always 2. This is enforced by the smart constructor(s).
newtype AABB = AABB (Matrix Double)

-- | A point in the physics universe is just a list of the object's
--   coordinates. The length of the vector is the number of dimensions its's
--   defined for.
--
--   Use the smart constructor 'point2D' to build points.
newtype Point = Point (Vector Double)

-- | "Oriented Bounding Box." A bounding box with any arbitrary orientation.
--
--   Use smart constructor 'obb2D' to build OBBs.
--
--   * Implementation Notes:
--
--   An oriented bounding box consists of a rotation/translation matrix,
--   and a vector of n-dimensional extents. For example, a 2x2 bounding box
--   would have the extents [1,1].
--
--   We take a simplified approach in 2D - we just use the four corners of the
--   box, specified in a clockwise orientation. This turns our point-OBB
--   intersection test into 8 vector subtractions, and 4 dot products. Needles
--   to say, that's ridiculously cheap.
data OBB = OBB2D {-# UNPACK #-} !Point
                 {-# UNPACK #-} !Point
                 {-# UNPACK #-} !Point
                 {-# UNPACK #-} !Point
         | OBB !(Matrix Double) -- The transformation matrix.
               !(Vector Double) -- The extent matrix.

-- | Builds a 2D point out of a tuple.
point2D :: Integral a => (a, a) -> Point
point2D (x, y) = Point . fromList $ map fromIntegral [x, y]

-- | Constructs a 2-dimensional axis-aligned bounding box.
aabb2D :: (Double, Double) -> (Double, Double) -> AABB
aabb2D (xa, xb) (ya, yb) = AABB $ (2 >< 2) [ min xa xb, max xa xb
                                           , min ya yb, max ya yb
                                           ]

-- | Constructs a 2D oriented bounding box from a list of verticies, specified
--   in a clockwise direction.
--
--   Seriously - clockwise. Don't fuck this up. It's _very_ important. No
--   checks are done.
--
--   Clockwise.
obb2DCorners :: Point -> Point -> Point -> Point -> OBB
obb2DCorners = OBB2D

-- | The type class for our 'intersects' method. This lets us perform double
--   dispatch on the types of both of intersects' arguments.
class Intersectable a b where
    -- | We define intersects for as many pairs of geometry as possibly,
    --   allowing us to do double-dispatch collision detection.
    intersects :: a -> b -> Bool

between :: Ord a => (a, a) -> a -> Bool
between (l, r) x = l < x && x < r
{-# INLINE between #-}

instance Intersectable AABB AABB where
    -- | Two AABBs intersect if any of their edges lie between the other's
    --   edges.
    intersects (AABB ba) (AABB bb) = any testDim . zip (toRows ba) $ toRows bb
        where
            testDim :: (Vector Double, Vector Double) -> Bool
            testDim (a, b) = let a'@(lx, hx) = (a @> 0, a @> 1)
                                 b'@(ly, hy) = (b @> 0, b @> 1)
                              in  between a' ly || between a' hy
                               || between b' lx || between b' hx
    {-# INLINE intersects #-}

instance Intersectable AABB Point where
    intersects (AABB bs) (Point p) = all inRange . zip (toRows bs) $ toList p
        where
            inRange :: (Vector Double, Double) -> Bool
            inRange (bnds, x) = between (bnds @> 0, bnds @> 1) x
    {-# INLINE intersects #-}

instance Intersectable Point AABB where
    intersects x y = intersects y x
    {-# INLINE intersects #-}

obbBadInput :: a
obbBadInput = error "Invalid dimensions in the OOB intersection vectors being tested."

instance Intersectable OBB Point where
    intersects (OBB2D (Point w) (Point x) (Point y) (Point z)) (Point p)
                   | any (/= 2) [ dp, dw, dx, dy, dz ] = obbBadInput
                   | otherwise = all (>= 0) [ (p - w) <.> (x - w)
                                           , (p - x) <.> (y - x)
                                           , (p - y) <.> (z - y)
                                           , (p - z) <.> (w - z)
                                           ]
        where
            dp = dim p
            dw = dim w
            dx = dim x
            dy = dim y
            dz = dim z
    intersects (OBB lm exts) (Point p) | allEqual [ dlmx, dlmy, dp ] = obbBadInput
                                       | otherwise = intersects aabb transformed
        where
            transformed :: Point
            transformed = Point $ inv lm <> p
            aabb :: AABB
            aabb = AABB . fromRows . map (\e -> fromList [-e, e]) $ toList exts

            (dlmx, dlmy) = (cols lm, rows lm)
            dp           = dim p

-- | Are all elements in the list equal?
allEqual :: Eq a => [a] -> Bool
allEqual []     = True
allEqual (x:xs) = fst $ foldl' (\(b, v) v' -> (b && v == v', v)) (True, x) xs
-- we use foldl', since in most cases, the elements _will_ be equal. If the
-- converse were true, I would prefer foldr since && is lazy.

instance Intersectable Point OBB where
    intersects x y = intersects y x
    {-# INLINE intersects #-}

max_stars_repo_name: bfops/Chess

max_issues_repo_name: bfops/Chess

max_forks_repo_name: bfops/Chess


text:
{-# LANGUAGE TupleSections ,ViewPatterns ,NoMonomorphismRestriction ,BangPatterns ,ScopedTypeVariables #-}
module Main
where
{- Test performance of encoding against other packages -}

{-# LANGUAGE MultiParamTypeClasses ,DeriveGeneric ,DeriveDataTypeable ,ScopedTypeVariables ,GADTs ,NoMonomorphismRestriction ,DeriveGeneric ,DefaultSignatures ,TemplateHaskell ,TypeFamilies ,FlexibleContexts ,CPP #-}

-- cabal update; cabal install binary binary-bits cereal buffer-builder bytes bytestring bits

{-
Compile:
cd /Users/titto/workspace/quid2;cabal install

Run plain, measure time externally:
cd /Users/titto/workspace/quid2;cabal install;time /Users/titto/workspace/quid2/.cabal-sandbox/bin/benchsmall

cd /Users/titto/workspace/quid2;cabal install;time /Users/titto/workspace/quid2/.cabal-sandbox/bin/benchsmall +RTS -sstderr -h;hp2ps -c benchsmall

Profiling
cd /Users/titto/workspace/quid2;cabal install --enable-executbitEncoder :: Int -> Encoding -> L.ByteString
bitEncoder numBits (Writer op) = L.fromStrict $ bitEncoderStrict numBits op
able-profiling -p;/Users/titto/workspace/quid2/.cabal-sandbox/bin/benchsmall +RTS -ph -sstderr;cat benchsmall.prof;hp2ps -c benchsmall;profiteur benchsmall.prof;open file:///Users/titto/workspace/quid2/benchsmall.prof.eUnsigned64 :: Word64 -> Encodinghtml

Benchmarking
cd /Users/titto/workspace/quid2;cabal install;/Users/titto/workspace/quid2/.cabal-sandbo# - store-0.4.1
# - store-core-0.4
x/bin/benchsmall

Benchmarking using cbor test suite, see LEGGIMI

Profiling using eventlog:
cd /Users/titto/workspace/quid2;cabal install;time /Users/titto/workspace/quid2/.cabal-sandbox/bin/benchs-- mall +RTS -l;ghc-event-- let v = [1::Word64 .. 4000000]
  -- evaluate $ force' (NF v)
  -- putStrLn . commas . length $ v

  s-analyze benchsmall.eventlog;cat benchsmall.totals.txt

-}

{-
Decoding
|               | ln2 |
| binary        | 23  |
| BitEncoderE   | 23  |

DecM
|                    | treeN μs         | ln2 μs | ln3 large ms | tree3 ms| tuple μs | One ns
| new-binary Custom  | 76               |
| binary 0.7.2.2     | 48               |  83    | 150         | 177
| cereal 0.4.1.1     | 58               |        |
| quid2 "binary-bits"| 87               |        |  271          | 426
| .. getBoolFast     | 78         wordsT = ("words",wordsV)
wordsV = (18::Word,33::Word8,1230::Word16,9990::Word32,1231232::Word64)
      |        |
| .. Custom        -ddump-opt-cmm  |       -ddump-opt-cmm           |        |
| .. DECODER_BITS    | 50               |        | 256:156-306   | 282:248-312
| .. "binary-strict" | 370              |        | 1052          | 1275
| .. Cstm BitDecoder
| DEC [1]            | 41               | 123    | 95            | 167
| DEC [2]            | [42]             | 88     | [60]          | 167
| DEC [3]            | 40               | [78]   | 132/60        | [163]
| .. DEC1            | 73               | 10     | 232:152-268   | 209:181-232
| .. runDec3         |                  |        | 256:153-300!  | 235:213-250
| .. runDec33        |                  |        | 261:194-304!  | 249:221-273
| .. runDec2         |                  |        | 266           | 238
| .. runDec22        |                  |        | 236:157-279   | 220:213-234
| .. DEC2 custom     |                  |        | 225           | 387
| .. .. Strict       | 56               | 12     |

[1] -DDEC=2 -D__INLINE -DDECODER_DATA -DGOCASES=0 -DDECODER_QUID2

[2] -DDEC=2 -D__INLINE -DDECODER_PRIM -DGOCASES=0 -DDECODER_QUID2

[3] -DDEC=2 -D__INLINE -DDECODER_GETTER -DGOCASES=0 -DDECODER_QUID2

EncDec
|                    | treeN μs| ln2 μs  | ln3 ms | IO ln3 ms | tuple μs | One ns
| new-binary Custom  | 91      | 240     | 25     |  15      | 1.7     | 388
| binary 0.7.2.2     | 57      | 164     | 16-17  |  7-8     | 1       | 562
| .. 0.7.2.3-0.7.4.0 |         |         | !105!  |
| cereal             | 70      | 190     | 16     |  8-9     | 1.2     | 493
| BitEncoderD        | 82      | 239     | 18     |  13-14   | 1.2     | 458
| BitEncoderE Custom | 74      | 201     | 18     |          | 1.2     | 439
| BitEncoderE        | 88      | 248     | 20     |          | 1.5     | 528
| BitEncoderF        |                   | 19

Enc
|                    | treeN | ln2   | ln3 ms | tuple    | One
| new-binary Custom  | 28    | 67    | 3.5    |  655     | 223
| binary 0.7.2.2     | 29    | 89-92 | 4.5    |  563-592 | 315
| cereal             | 28    | 81    | 4      |  551     | 335
| BitEncoderD        | [26]  | 69-75 | 3.8    |  357     | 224
| BitEncoderD Custom | 22-24 | 63-64 | 3.1    |  357     | 224
| BitEncoderE        | 30    | 82    | 4.3    |  585     | 242
| !BitEncoderF       | 27    | [65]  | [3.2]  |  444     | 220
| .. Custom          | 24-25 | [51]  | [2.5]  |  371     | 224
| .. Custom Strict   | 21    | 63    | [10.19]|  353     | 218
| ..Cst toBitBuilderE| 80    | [54]  | [2.86] |  431     | 214

|           | tstEnc lBool | tstEnc lN2 | tstEncodeDecode lBool | tstEncodeDecode lN2 | tstDec lN2
new-binary  | 6.8          | 65-69      | 22.7                  | 232-237             |
binary      |              | 86-88      | 13.7                  | 161-169
binary custom |            | 83-85
cereal      |              | 80-82          | 15.7              | 185-190
ByteBuilder |              | 78            | 20.8               | 236
Data.Binary.Quid2 BitEnc   | 135-141             |
.. BitEncoderF Custom       | 55-56            |
.. BitEncoderE              | 59-60                             | 230                  | 165
.. BitEncoderE Custom       | 51-53        |                    | 160-195              | 136
.. ByteEncoderF             | 121-135      |
.. ByteEncoderF Custom      | 55           |
Byte        |              | 158           | 35                  | 397-413
Byte2       |              | 187           | 37.5                | 412
Bit2 (Coding 64) | 20.8    |               |
Bit2 (Coding 8)  | 22      |               |
Bit  (BitPut orig) | 15-17 | 384           | 55                   | 831-851
Bit3 (BitPut 8)    | 14-17 | 220                                  | 664-690
Bit3 (BitPut 64)   | 15-18 | 242
-}

-- ghc -O2 --make Bench;./Bench
import           Control.DeepSeq
import Criterion.Main
import Criterion.IO
import Criterion.Types
import Types
import Test.Data
import           Test.Data.Values
import PkgBinary
import qualified PkgBinary as B
-- import PkgCereal
import qualified PkgCereal as CC
import PkgCBOR
import qualified PkgCBOR as C
import PkgFlat
import qualified PkgFlat as F
import PkgStore
import qualified PkgStore as S
import qualified Data.Binary.Serialise.CBOR as CBOR
import qualified Data.ByteString.Lazy as B
import qualified Data.ByteString.Lazy as L
import qualified Data.ByteString as BS
import qualified Data.Text as T
import Control.Exception
import Control.Monad
import Control.Applicative
--import System.Directory
import Debug.Trace -- (traceEventIO)
import           Criterion.IO
import           Criterion.Types
-- import           Statistics.Resampling.Bootstrap
import qualified Data.Map as M
import qualified Data.Flat as F
import qualified Data.Flat.Encoder as F
import System.Directory
import Report
import Weigh

event label =
  bracket_ (traceEventIO $ "START " ++ label)
           (traceEventIO $ "STOP "  ++ label)

t = mainBench

-- tl1 = [tstLen PkgBinary tree1
--       ,tstLen PkgBinary tree2
--       ,tstLen PkgBinary lBool
--       --,tstLen PkgBinary lN
--       ]

-- tl2 = [--tstLen PkgFlat tree1
--       -- ,tstLen PkgFlat tree2
--       tstLen PkgFlat lBool
--       --,tstLen PkgFlat lN
--       ]

l = [tstSize tupleT,tstSize tupleBools,tstSize tupleWords,tstSize arr0, tstSize arr1, tstSize asciiStrT,tstSize unicodeStrT,tstSize unicodeTextT , tstSize lN3T, tstSize treeNLargeT]
-- [("tuple",(8.0,3.0,3.0)),("tupleBools",(7.0,4.5,4.5)),("tupleWord",(1.5555555555555556,3.5555555555555554,3.5555555555555554)),("[Bool]",(3.999992000032,4.000015999936,4.000015999936)),("[Word]",(1.3416352722052283,4.691580498057688,4.691580498057688)),("asciiStr",(0.8888925432066282,0.8888952098709245,3.555559506169328)),("unicodeStr",(0.9322290145814176,0.932230773922345,2.345792594816747)),("unicodeText",(0.9960948808600285,0.9960967607312449,1.2532524905160498)),("lN3",(9.41174726647702,4.70587363323851,4.70587363323851)),("tree3",(10.909087272727273,5.454543636363637,5.454543636363637))]
-- tstSize :: (CBOR.Serialise a, Flat a) => (String, a) -> (String,Double)
tstSize nv@(n,v) = (n,) $ (tstLen PkgCBOR v,tstLen PkgBinary v,tstLen PkgStore v)

tstLen:: (Serialize c a, Flat a1) => (a1 -> c a) -> a1 -> Double
tstLen k v = let l0 = lser PkgFlat v
                 l1 = lser k v
             in (l1 / l0)

-- lser :: (Serialize c1 a1) => (a -> c1 a1) -> a -> Double
lser k = fromIntegral . Prelude.length . ser . k

-- main = mainWeight
main = mainB

{-# RULES "cl1" forall a. cl [a] = 1 ; #-}
{-# RULES "cl2" forall a b. cl [a,b] = 2 ; #-}

{-# INLINE [0] cl #-}
-- {-# NOINLINE cl #-}
cl :: [a] -> Int
cl l = error "BAD" -- length l

yv = cl [2,4]
xv = yv + cl [2]

-- main = mainEvents
mainB = do
  forceCafs

  prints
  mainBench

  --mainProfile
  -- mainAlternatives

prints = do
  -- plf treeNLargeT
  -- -- plf treeNLargeT >> plf treeNNNLargeT >>
  -- plf nativeListT

  -- print (yv,xv)
  -- print $ pt carT

  print $ F.getSize Three
  print $ F.getSize (Leaf Three)
  print $ F.getSize (treeNLarge)
  print $ F.getSize (treeNNNLarge)
  print $ F.getSize (nativeList)

  -- print $ B.unpack (F.flat v1)
  -- print $ B.unpack (F.flat v2)
  -- print $ B.unpack (F.flat t33)
  -- print $ [33] == B.unpack (F.flat v1)
  -- print $ [89] == B.unpack (F.flat v2)
  -- print $ [220,221,25] == B.unpack (F.flat t33)

  print $ 1175000 == B.length (F.flat treeNNNLarge)
  print $ 550000 == B.length (F.flat treeNLarge)
  print $ 2125001== B.length (F.flat nativeList)
  print $ (B.length (F.flat unicodeTextT),B.length (F.flat unicodeTextUTF8T),B.length (F.flat unicodeStrT))
  pp car1
  -- ),B.length (serialize PkgStore car1))
  -- print $ 1 == B.length (F.flat v2)
  -- print $ 3 == B.length (F.flat t33)
  -- print $ 10 == B.length (F.flat wordsV)
  -- print $ 67 == B.length (F.flat car1)
  -- --print F.eee
  print $ ("flat",encodeOnly PkgFlat car1,"store",encodeOnly PkgStore car1,"cbor",encodeOnly PkgCBOR car1)

  pp (Node (Leaf One) (Leaf Two))
  pp ((One,Two,Three),(One,Two,Three),(One,Two,Three))
  pp t33
  pp (asN33 4)
  pp v1
  pp v2
  pp [Five,Five,Five]
  pp "aaa"
  pp vw
  pp wordsV
  --nprint (F.encode (undefined::Various))
    where pp v = print (v,F.encode v) --  F.<> F.eFiller)

mainWeight = do
  forceCafs
  mainWith $ do
    --ww nativeListT
    ww treeNLargeT
    --ww wordsT

  pp nativeListT
  pp treeNLargeT
  pp wordsT
  pp unicodeStrT
  --pp unicodeTextUTF8T

  where
    ww o = do
      weightT "store" PkgStore o
      weightT "binary" PkgBinary o
      -- weightT "cbor" PkgCBOR o
      weightT "flat" PkgFlat o

pp o = do
  pl "store" PkgStore o
  pl "binary" PkgBinary o
  pl "flat" PkgFlat o

plf = pl "flat" PkgFlat

pl kname k (n,v) = putStrLn (unwords ["Len serialise",kname,n,commas $ encodeOnly k v])
weightT kname k (n,v) = func (unwords ["serialise",kname,n]) (encodeOnly k) v

encodeOnly k = B.length . serialize . k

mainEvents = mainAlternatives

-- #define PKG_BINARY

mainProfile = mapM_ (\_ -> enc treeNT) [1..1000] -- mapM_ (\_ -> enc (0x34::Word) >> enc tree1 >> enc tree2) [1..10]
  where enc = evaluate . B.length . serialize . PkgFlat -- PkgBit2

mainAlternatives = dec_ tupleT >> dec_ tupleBools >> dec_ arr0 >> dec_ arr1 >> dec_ asciiStrT >> dec_ unicodeStrT >> dec_ lN3T >> dec_ treeNLargeT -- dec_ sbs >> dec_ lbs >>
  where
    -- dec1 = evaluate . and . Prelude.map (\_ -> des) $ [1..1]
    -- des = let (Right (PkgFlat a)) = deserialize PkgFlat.serlN2 in (a == lN2)
    dec_ dt@(n,_) = event (Prelude.concat [testName,"_",n]) $ (mapM (\_ -> d dt) [1..1])
#ifdef PKG_BINARY
    d = op "binary" PkgBinary
#elif PKG_QUID2
    d = op "quid2" PkgFlat
#elif PKG_STORE
    d = op "store" PkgStore
#elif PKG_CBOR
    d = op "cbor" PkgCBOR
#endif
#if OP_DEC
    testName = "dec"
    op = decM
#elif OP_ENC
    testName = "enc"
    op _ pkg v = (evaluate . B.length . serialize . pkg $ v)
#elif OP_ENCDEC
    testName = "encDec"
    op = encDecM
#endif

workDir = "."

mainBench = do
  mainBench_ (reportsFile workDir)
  ms <- updateMeasures_ workDir

  printMeasuresDiff ms
  -- printMeasuresAll ms
  printMeasuresCurrent ms

--mainBench = defaultMainWith (defaultConfigFilePath {
mainBench_ jsonReportFile = defaultMainWith (defaultConfig {jsonFile= Just jsonReportFile}) (
  concat [
   --tstDec carT
   --tstDec nativeListT
   --,tstDec treeNLargeT
   --,tstDec treeNNNLargeT
   --,tstDec wordsT,tstDec words0T
   --,tstDec vwT,tstDec vfT,tstDec viT
   -- ,tstDec v2T
   -- -- ,tstDec charT
   -- ,tstDec unicharT
   -- ,tstDec lN2T
   -- -- --,tstDec lN3T
   -- -- flat fails to complete:
   --tstDec treeN33LargeT
   tstDec mapT--,tstDec seqNT
   -- ,tstDec asciiStrT
   --,tstDec unicodeStrT
   --,tstDec unicodeTextT
   --,tstDec unicodeTextUTF8T
   -- --,tstDec sbs,tstDec lbs
   --,tstDec shortbs
   ]

   -- ++ [
   --    tstMaxSize treeNLargeT
   --    ,tstMaxSize treeNNNLargeT,tstMaxSize nativeListT
   --    ,tstMaxSize wordsT,tstMaxSize words0T
   --    ,tstMaxSize vwT,tstMaxSize vfT,tstMaxSize viT
   --    ,tstMaxSize carT,tstMaxSize seqNT
   -- ]
     -- -- ,tstEnc viiT,,tstEnc intsT
  -- ,tstEnc sbs,tstEnc lbs --,tstEnc shortbs

  --,tstEnc treeN33LargeT
  --,tstEnc floatT,tstEnc doubleT
  --   ++ join [tstDec carT
  --           --,tstDec wordsT,tstDec vwT
  --           --,tstDec lN3T,tstDec treeNLargeT,tstDec lN2T
  --
  --]
  )

{-
   tstEnc (0x34::Word)
  ,tstEnc (0x1289::Word)
  ,tstEnc (0x221289::Word)
  ,tstEnc (0x55221289::Word)
  ,tstEnc (0x1233456712345533::Word)

  ,tstEnc ((-120)::Int)
  ,tstEnc ((-1203321312)::Int)

   tstEnc (3.4512312::Float)
  ,tstEnc (-0.4511231232312::Double)

  ,tstEnc 'z'
  ,tstEnc '中'

  -- sequences
  ,tstEnc hl1
  ,tstEnc hl2
  ,tstEnc hl3

  ,tstEnc s1
  ,tstEnc s2
  ,tstEnc s3
  ,tstEnc s4

  ,tstEnc t1
  ,tstEnc t2
  ,tstEnc t3
  ,tstEnc t4

  ,tstEnc b1
  ,tstEnc b2
  -- ,tstEnc lb1
  -- ,tstEnc lb2

  -- generic
  -- ,tstEnc True
  -- ,tstEnc ('a','b',('c','d',('e','f','g')))
-}
  --

  -- tstEnc treeNT,tstEnc treeNLargeT,tstEnc lN2T,tstEnc lN3T

  -- ,tstEnc oneT

  -- recursive structure  tstEnc lBool

  -- ,tstEnc tree2

  --
  -- ,tstDecM lN2T
  -- tstEncDecM lN2T,tstEncDecM lN3T
  -- tstEncDecM treeNT,tstEncDecM treeNLargeT,tstEncDecM lN2T,tstEncDecM lN3T,
  -- tstDecM treeNT,tstDecM lN2T,tstDecM lN3T,tstDecM treeNLargeT

  -- tstEncDecM tuple0T,tstEncDecM tupleT,
  -- tstEncDecM word8T,tstEncDecM int8T
  --tstEncDecM word64T,tstEncDecM int64T,tstEncDecM integerT,

  --tstEnc carT
  -- ,tstEnc floatT,tstEnc doubleT,tstEnc floatsT,tstEnc floatsUnaT,tstEnc vfT

  -- tstEnc sbs,tstEnc lbs
  --,tstEnc shortbs
  --,tstEnc wordsT,tstEnc vwT
  --,tstEnc intsT,tstEnc viT,tstEnc viiT
  --tstEnc charT,tstEnc unicharT
  -- ,tstEnc nativeListT,tstEnc unicodeStrT,tstEnc treeNLargeT
   -- ,tstEnc lN3T
  --,tstEnc asciiStrT

  -- ,tstEnc unicodeTextT

  -- ,tstEnc treeN33LargeT
  --,tstEnc lN3T

    -- tstEnc v2T,
 -- tstEnc t33T,
 -- tstEnc treeVariousT
 -- ,tstEnc treeNT

  --tstEnc carT,tstEnc floatsT,tstEnc tupleT,tstEnc tupleBools,tstEnc tupleWords,tstEnc arr0, tstEnc arr1, tstEnc asciiStrT,tstEnc unicodeStrT,tstEnc unicodeTextT , tstEnc lN3T,
  --tstEncDec floatT,tstEncDec doubleT,tstEncDec tupleT,tstEncDec tupleBools,tstEncDec tupleWords,tstEncDec arr0, tstEncDec arr1, tstEncDec asciiStrT,tstEncDec unicodeStrT,tstEncDec unicodeTextT , tstEncDec lN3T, tstEncDec treeNLargeT
  -- tstEncDecM tupleT,tstEncDecM tupleBools,tstEncDecM tupleWords,tstEncDecM arr0, tstEncDecM arr1, tstEncDecM asciiStrT,tstEncDecM unicodeStrT,tstEncDecM unicodeTextT , tstEncDecM lN3T, tstEncDecM treeNLargeT
  --tstEncDecM treeNT,tstEncDecM treeNLargeT,tstEncDecM lN2T,tstEncDecM lN3T
  -- ,tstDecM treeNT,tstDecM treeNLargeT,tstDecM lN2T,tstDecM lN3T
  -- ,tstDecM treeNLargeT

  --tstED lN3T,tstED treeNLargeT

  --,tstEncDecM treeNT
{-
  ,tstEncDec treeN
  ,tstEncDec lN2
  ,tstEncDec lN3
  ,tstEncDec (Two,One,(Five,Three,(Three,(),Two)))
  ,tstEncDec One
-}
 -- ,tstEncodeDecode "BEPkg" BEPkg,  tstEncodeDecode "PkgCBOR" PkgCBOR
{-
  ,tstEncodeDecode "PkgBinary" PkgBinary
  -- ,tstEncodeDecode "PkgCereal" PkgCereal
  ,tstEncodeDecode "PkgFlat" PkgFlat
  -- ,tstEncodeDecode "PkgByte" PkgByte
  -- ,tstEncodeDecode "PkgByte2" PkgByte2
  -- ,tstEncodeDecode "PkgBit" PkgBit
  -- ,tstEncodeDecode "PkgBit2" PkgBit2
-}
{-
  ,tstEncodeDeep "BEPkg" BEPkg BEPkg
  ,tstEncodeDeep "PkgBinary" PkgBinary PkgBinary
  ,tstEq
-}
                                                               --(join [
tstED t = bgroup ("tstEncDecsM") [tstEnc t
                                 --,tstEncM t
                                 ,tstDecM t
                                 ,tstEncDecM t,tstEncDec t]

tstEncDecM lv@(vn,v) =
  bgroup ("tstEncDecM") [
              tst "quid2" PkgFlat
             ,tst "cbor" PkgCBOR
             -- ,tst "new-binary" BEPkg
             ,tst "binary" PkgBinary
             ,tst "store" PkgStore
             --,tst "cereal" PkgCereal
             ]
      where
        tst n c = bench (unwords[n,vn]) $ nfIO (encDecM n c lv)

df n vn = Prelude.concat[workDir,"/benchmarks/data/tstEncDec-",vn,"-",n]

encDecM n c (vn,v) = do
          let f = df n vn
          B.writeFile f . serialize . c $ v
          rs <- deserialize . B.fromStrict <$> BS.readFile f
          -- print rs
          let Right r = rs
          when (r /= c v) $ dataErr f
          -- r <- BS.readFile f
          -- evaluate (BS.length r)
          -- evaluate (B.length . serialize . c $ v)

decM n c (vn,v) = do
          let f = df n vn
          Right r <- deserialize . B.fromStrict <$> BS.readFile f
          when (r /= c v) $ dataErr f

dataErr f = error . unwords $ ["Data mismatch on",f]

tstEncM (vn,v) =
  bgroup ("tstEncM") [
              tst "quid2" PkgFlat
             -- tst "new-binary" BEPkg
             ,tst "binary" PkgBinary
             --,tst "cereal" PkgCereal
             ]
      where
        tst n c = bench (unwords[n,vn]) $ nfIO (rw n c)
        rw n c = do
          let f = df n vn
          removeFile f
          B.writeFile f . serialize . c $ v
          -- Right r <- deserialize . B.fromChunks . (:[]) <$> BS.readFile f
          -- when (r /= c v) $ print "Data mismatch"
          r <- BS.readFile f
          evaluate (BS.length r)
          -- print $ BS.length r
          -- evaluate (B.length . serialize . c $ v)

tstDecM d@(vn,v) =
  bgroup ("tstDecM") [
              tst "quid2" PkgFlat
             -- ,tst "new-binary" BEPkg
             ,tst "cbor" PkgCBOR
             ,tst "binary" PkgBinary
--             ,tst "binary-strict" PkgBinaryStrict
             --,tst "cereal" PkgCereal
             ]
      where
        tst n c = bench (unwords[n,vn]) $ nfIO (decM n c d)
{-
        rw n c = do
          let f =
          Right r <- deserialize . B.fromStrict <$> BS.readFile f
          when (r /= c v) $ print "Data mismatch"
          -- r <- BS.readFile f
          -- when (BS.length r < 0) $ print "Data mismatch"
-}

tstEncDec (vn,v) = let nm s = unwords[s,vn] -- unwords [s,Prelude.take 400 $ show v]
 in bgroup ("tstEncDec") [
  bench (nm "flat")   $ nf (encodeInverseOfDecode PkgFlat) v
  ,bench (nm "cbor") $ nf (encodeInverseOfDecode PkgCBOR) v
  ,bench (nm "binary") $ nf (encodeInverseOfDecode PkgBinary) v
  ,bench (nm "store") $ nf (encodeInverseOfDecode PkgStore) v
  --,bench (nm "cereal")     $ nf (encodeInverseOfDecode PkgCereal) v
  ]

encodeInverseOfDecode k v = let x = k v
                                Right r = (deserialize . serialize) x
                            in x == r

tstEncE (vn,v) =
  let nm s = Prelude.concat [vn,"-",s] -- unwords [s,Prelude.take 400 $ show v]
  in bgroup ("tstEnc") [
    --bench (nm "store")     $ nf (encodeOnly PkgStore) v
    --,bench (nm "binary")     $ nf (encodeOnly PkgBinary) v
    --,bench (nm "cbor")     $ nf (encodeOnly PkgCBOR) v
    bench (nm "flat")  $ nf (encodeOnly PkgFlat) v
    ]
  where encodeOnly k = B.length . serialize . k

-- tstDec :: forall b. (S.Store b, F.Flat b,Eq b) => (String, b) -> Benchmark

tstMaxSize (vn,v) =
  let nm s = Prelude.concat [vn,"-",s] -- unwords [s,Prelude.take 400 $ show v]
  in bgroup ("tstMaxSize") [
    bench (nm "flatClassSizeSimple")  $ nf (F.getSize) v
    --,bench (nm "store")  $ nf (S.getSize) v
    ]
  where encodeOnly k = B.length . serialize . k

tstEnc (vname,obj) =
  let nm s = Prelude.concat [vname,"-",s]
  in bgroup ("tstEnc") $ map (\(pkg,_,s,d) -> bench (nm pkg) $ nf (L.length . s) obj) pkgs

tstDec (vname,obj) =
  let nm s = Prelude.concat [vname,"-",s] -- unwords [s,Prelude.take 400 $ show v]
  in [
     bgroup ("tstEncDec") $ map (\(_,pkg,s,d) -> bench (nm pkg) $ nf (\obj -> d (s obj) == Right obj) obj) pkgs
     ,bgroup ("tstEnc") $ map (\(pkg,_,s,d) -> bench (nm pkg) $ nf (B.length . s) obj) pkgs

     --bgroup ("tstDec") $ map (\(_,pkg,s,d) -> env (return $! s obj) $ (\bs -> bench (nm pkg) $ nfIO ((Right obj ==) <$> (return . force . d) bs))) pkgs
     -- in [bgroup ("tstEncDec") $ map (\(pkg,s,d) -> bench (nm pkg) $ nf (\obj -> d (s obj) == Right obj) obj) pkgs
     --,bgroup ("tstDec") $ map (\(pkg,s,d) -> env (return $ s obj) $ (\bs -> bench (nm pkg) $ nf (\obj -> Right obj == d bs) obj)) pkgs
     --,bgroup ("tstDec") $ map (\(pkg,s,d) -> env (return $ s obj) $ (\bs -> bench (nm pkg) $ nfIO ((Right obj ==) <$> return (d bs)))) pkgs
      --bgroup ("tstDec") $ map (\(_,pkg,s,d) -> env (return $ s obj) $ (\bs -> bench (nm pkg) $ whnfIO ((Right obj ==) <$>  (return . force . d $ bs)))) pkgs
     --,bench ("tstEq/"++nm "any") $ nf (\o -> o == o) obj
     ]

pt (n,v) = map (\(_,pkg,s,d) -> Right v == d (s v)) pkgs

--pkgs :: (CC.Serialize a ,C.Serialise a,S.Store a,B.Binary a,F.Flat a) => [(String,String,a -> L.ByteString,L.ByteString -> Either String a)]
-- pkgs = [S.sd,B.sd,C.sd,CC.sd,F.sd]

-- pkgs :: (C.Serialise a) => [(String,String,a -> L.ByteString,L.ByteString -> Either String a)]
-- pkgs = [C.sd]

-- pkgs :: (C.Serialise a,S.Store a,B.Binary a,F.Flat a) => [(String,String,a -> L.ByteString,L.ByteString -> Either String a)]
-- pkgs = [S.sd,B.sd,C.sd,F.sd]

pkgs :: (B.Binary a,S.Store a,F.Flat a) => [(String,String,a -> L.ByteString,L.ByteString -> Either String a)]
pkgs = [S.sd,F.sd]

-- pkgs :: (F.Flat a) => [(String,String,a -> L.ByteString,L.ByteString -> Either String a)]
-- pkgs = [F.sd]

tstEncodeDeep name k1 k2 = bgroup ("serialize " ++ name) [
   bench "tree1"  $ nf (encodeOnly k1) tree1
  ,bench "tree2" $ nf (encodeOnly k2) tree2
  ]
  where encodeOnly k = B.length . serialize . k

tstDec2 = bgroup ("tstDec2") [
  bench "lN2-binary" $ nf (\v -> let Right (PkgBinary a) = deserialize PkgBinary.serlN2 in a == v) lN2
  ,bench "lN2-flat" $ nf (\v -> let Right (PkgFlat a) = deserialize PkgFlat.serlN2 in a == v) lN2
  ]

tstEncodeDecode name k1 = bgroup ("serialize+deserialise " ++ name) [
  -- bench "tree1" $ whnf (encodeInverseOfDecode k1) tree1
  --,bench "tree2" $ whnf (encodeInverseOfDecode k2) tree2
   -- ,bench "car1" $ whnf encodeInverseOfDecode car1
   bench "lN2" $ nf (encodeInverseOfDecode k1) lN2
  ,bench "lN3" $ nf (encodeInverseOfDecode k1) lN3
  -- bench "treeN" $ nf (encodeInverseOfDecode k1) treeN
  ]

-- dec p b v = let (Right p) = deserialize b in (unPkg p == v)

-- encodeInverseOfDecode :: (Serialize a, Eq a) => (t -> a) -> t -> Bool
{-
tstEq = bgroup ("equal") [
  bench "lBool" $ whnf eq lBool
  ,bench "tree1" $ whnf eq tree1
  ,bench "tree2" $ whnf eq tree2
  ,bench "car1" $ whnf eq car1
  ]
-}

eq v = (v == v) == True

max_stars_repo_name: AnisimoffNikita/flat

max_issues_repo_name: AnisimoffNikita/flat

max_forks_repo_name: AnisimoffNikita/flat


text:
{-# LANGUAGE FlexibleContexts    #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE PolyKinds           #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
-- |
-- Module      : Bench.Util
-- Copyright   : [2017..2020] Trevor L. McDonell
-- License     : BSD3
--
-- Maintainer  : Trevor L. McDonell <trevor.mcdonell@gmail.com>
-- Stability   : experimental
-- Portability : non-portable (GHC extensions)
--

module Bench.Util where

import Data.Complex
import System.Random.MWC

data ArgType (a :: *) = AT

showType :: forall proxy a. Show (ArgType a) => proxy a -> String
showType _ = show (AT :: ArgType a)

instance Show (ArgType Float)            where show _ = "Float"
instance Show (ArgType Double)           where show _ = "Double"
instance Show (ArgType (Complex Float))  where show _ = "Complex Float"
instance Show (ArgType (Complex Double)) where show _ = "Complex Double"

instance Variate e => Variate (Complex e) where
  uniform    gen = (:+) <$> uniform gen <*> uniform gen
  uniformR r gen =
    let (ur:+ui,vr:+vi) = r
    in  (:+) <$> uniformR (ur,vr) gen <*> uniformR (ui,vi) gen

infixr 0 $$
($$) :: (b -> a) -> (c -> d -> b) -> c -> d -> a
(f $$ g) x y = f (g x y)

max_stars_repo_name: statusfailed/accelerate-blas

max_issues_repo_name: statusfailed/accelerate-blas

max_forks_repo_name: statusfailed/accelerate-blas


text:
module Fourier.FFTSpec (main,spec) where

import Test.Hspec
import Test.QuickCheck()
import Data.Complex(Complex((:+)))
import Fourier.FFT(fft)

main :: IO ()
main = hspec spec

spec :: Spec
spec = do

  describe "FFT" $ do

    it "transform steady function" $ do
      let x = fft [0.0, 0.0, 0.0]
        in x `shouldBe` ([0.0, 0.0, 0.0])

    it "transform impulse function" $ do
      let x = fft [0.0, 0.0, 1.0, 0.0, 0.0]
        in x `shouldBe` ([
          1.0 :+ (-0.0),
          (-0.8090169943749473) :+ (-0.5877852522924732),
          0.30901699437494723 :+ 0.9510565162951536,
          0.30901699437494773 :+ (-0.9510565162951534),
          (-0.8090169943749477) :+ 0.5877852522924728])

    it "transform sinusoid function" $ do
      let x = fft [sin x | x <- [0,0.3 .. pi]]
        in x `shouldBe` ([
          6.654043728627056 :+ (-0.0),
          (-2.495308482972513) :+ (-0.3737795890019122)
          ,(-0.46503997943364506) :+ (-0.1328809251599644),
          (-0.1886713255075324) :+ (-7.112462151227344e-2),
          (-0.10350867238552199) :+ (-3.701914259223908e-2),
          (-7.449340401431373e-2) :+ (-1.1604714311766413e-2),
          (-7.449340401431312e-2) :+ 1.1604714311764401e-2,
          (-0.10350867238552014) :+ 3.701914259223764e-2,
          (-0.18867132550753074) :+ 7.112462151227178e-2,
          (-0.46503997943364067) :+ 0.13288092515996228,
          (-2.4953084829725105) :+ 0.3737795890019047])


text:
{-# LANGUAGE FlexibleContexts, NamedFieldPuns #-}

module School.Train.SimpleDescentUpdate
( simpleDescentUpdate ) where

import Numeric.LinearAlgebra (Container, Vector, add, scale)
import School.Train.TrainState (TrainState(..))
import School.Train.UpdateParams (UpdateParams)
import School.Types.PingPong (pingPongList, toPingPong)
import School.Unit.UnitParams (UnitParams(..))

update :: (Container Vector a, Num a)
       => a
       -> UnitParams a
       -> UnitParams a
       -> UnitParams a
update rate
       AffineParams { affineBias = bias
                    , affineWeights = weights}
       AffineParams { affineBias = biasGrad
                    , affineWeights = weightGrad } =
  AffineParams { affineBias, affineWeights } where
    affineBias = add bias (scale (-rate) biasGrad)
    affineWeights = add weights (scale (-rate) weightGrad)
update _ _ _ = EmptyParams

simpleDescentUpdate :: (Container Vector a, Num a)
                    => UpdateParams a
simpleDescentUpdate state@TrainState { paramDerivs
                                     , paramList
                                     , learningRate
                                     } = do
  newParamList <- toPingPong $
         zipWith (update learningRate)
                 (pingPongList paramList)
                 paramDerivs
  return state { paramList = newParamList }