text:
module STCR2Z1T0Edge where

import           Data.Array.Repa         as R
import           Data.Binary             (decodeFile)
import           Data.Complex
import           Data.List               as L
import           DFT.Plan
import           Filter.Utils
import           FokkerPlanck.MonteCarlo
import           FokkerPlanck.Pinwheel
import           Image.Edge
import           Image.IO
import           STC
import           System.Directory
import           System.Environment
import           System.FilePath
import           Types

main = do
  args@(numPointStr:numOrientationStr:sigmaStr:taoStr:lenStr:initialScaleStr:numTrailStr:maxTrailStr:theta0FreqsStr:thetaFreqsStr:histFilePath:pinwheelFlagStr:numIterationStr:writeSourceFlagStr:edgeFilePath:numThreadStr:_) <-
    getArgs
  print args
  let numPoint = read numPointStr :: Int
      numOrientation = read numOrientationStr :: Int
      sigma = read sigmaStr :: Double
      tao = read taoStr :: Double
      len = read lenStr :: Int
      initialScale = read initialScaleStr :: Double
      numTrail = read numTrailStr :: Int
      maxTrail = read maxTrailStr :: Int
      theta0Freq = read theta0FreqsStr :: Double
      theta0Freqs = [-theta0Freq .. theta0Freq]
      thetaFreq = read thetaFreqsStr :: Double
      thetaFreqs = [-thetaFreq .. thetaFreq]
      pinwheelFlag = read pinwheelFlagStr :: Bool
      numIteration = read numIterationStr :: Int
      writeSourceFlag = read writeSourceFlagStr :: Bool
      numThread = read numThreadStr :: Int
      folderPath = "output/test/STCR2Z1T0Edge"
  createDirectoryIfMissing True folderPath
  flag <- doesFileExist histFilePath
  radialArr <-
    if flag
      then R.map magnitude . getNormalizedHistogramArr <$>
           decodeFile histFilePath
      else do
        putStrLn "Couldn't find a Green's function data. Start simulation..."
        solveMonteCarloR2Z1T0Radial
          numThread
          numTrail
          maxTrail
          numPoint
          numPoint
          sigma
          tao
          1
          theta0Freqs
          thetaFreqs
          histFilePath
          (emptyHistogram
             [ (round . sqrt . fromIntegral $ 2 * (div numPoint 2) ^ 2)
             , L.length theta0Freqs
             , L.length thetaFreqs
             ]
             0)
  arrR2Z1T0 <-
    computeUnboxedP $
    computeR2Z1T0ArrayRadial
      radialArr
      numPoint
      numPoint
      1.5
      thetaFreqs
      theta0Freqs
  plan <- makeR2Z1T0Plan emptyPlan arrR2Z1T0
  xs <- parseEdgeFile edgeFilePath
  randomPonintSet <- generateRandomPointSet 40 numPoint numPoint
  let ys =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) -> (R2S1RPPoint (a - 150, b - 150, c, d)))
          xs
      zs =
        L.map
          (\(R2S1RPPoint (a, b, c, d)) ->
             (R2S1RPPoint (a - center numPoint, b - center numPoint, c, d)))
          randomPonintSet
      points = ys L.++ zs
  let bias = computeBiasR2T0 numPoint numPoint theta0Freqs points
      eigenVec =
        computeInitialEigenVectorR2T0
          numPoint
          numPoint
          theta0Freqs
          thetaFreqs
          points
  powerMethod1
    plan
    folderPath
    numPoint
    numPoint
    numOrientation
    thetaFreqs
    theta0Freqs
    arrR2Z1T0
    numIteration
    writeSourceFlag
    ""
    0.5
    bias
    eigenVec

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
import Statistics.LinearRegression
import Statistics.Sample
import qualified Data.Vector.Unboxed as U
import Control.Monad.Random
import Control.Monad
import Control.Applicative
import System.Random.MWC
import System.Random.MWC.Distributions
import qualified Data.Packed.Vector as V
import Graphics.Rendering.Plot

main = do
    mapM_ test [1..10]
    test_convergence
    test_robust
    test_variances

test k = do
    let n = 10000000
    let a = k*n + 1
    let b = (k+1)*n
    let xs = U.fromList [a..b]
    let ys = U.map (\x -> x*100 + 2000) xs
    putStrLn "linearRegression:"
    putStrLn . show $ linearRegression xs ys
    putStrLn "linearRegressionTLS:"
    putStrLn . show $ linearRegressionTLS xs ys

test_convergence = do
    let xs = U.fromList [1..10]
    let ys = U.fromList $ [1..5] ++ [0,0] ++ [8..10]
    let iter1 = linearRegression xs ys
    let ep = defaultEstimationParameters
    putStrLn "Initial iteration:"
    putStrLn . show $ iter1
    putStrLn "Successive iteration:"
    putStrLn . show $ converge ep xs ys iter1

getNormals :: Double -> Double -> Int -> IO [Double]
getNormals mean std n = do
    withSystemRandom . asGenIO $ \gen -> replicateM n (normal mean std gen)

testFigure :: U.Vector Double -> U.Vector Double -> (EstimatedRelation, EstimatedRelation, EstimatedRelation) -> Figure ()
testFigure xs ys (simple, non_robust, robust) = do
    let vxs = V.fromList . U.toList $ xs
    let vys = V.fromList . U.toList $ ys
    let dataset = (vxs, [   point vys Cross, line_func simple, line_func non_robust, line_func robust ])

    withTitle . setText $ "linreg test"
    setPlots 1 1
    withPlot (1,1) $ do
        addAxis XAxis (Side Lower) $ do
            setTicks Minor (TickNumber 5)
            withAxisLine $ do
                setLineWidth 1.0
        addAxis YAxis (Side Lower) $ do
            setTicks Minor (TickNumber 5)
            withAxisLine $ do
                setLineWidth 1.0
        setDataset dataset
        setRangeFromData XAxis Lower Linear
        setRangeFromData YAxis Lower Linear
        setLegend True NorthEast Inside
    where
        line_func (alpha,beta) = line ((\x -> alpha + beta*x) :: Function) (1.0 :: LineWidth)

test_robust = do
    putStrLn "generating random dataset for robust fit:"
    first_xs <- getNormals 0.0 10.0 800
    first_ys_errs <- getNormals  0.0 1.0 800
    let first_ys = zipWith (+) first_xs first_ys_errs
    last_xs <- getNormals  50.0 (sqrt 50) 200
    last_ys <- getNormals  0.0 (sqrt 50) 200
    let xs = U.fromList $ first_xs ++ last_xs
    let ys = U.fromList $ first_ys ++ last_ys
    putStrLn "robustFit test results:"
    (simple,non_robust,robust) <- evalRandIO (randTest xs ys)
    putStrLn "linearRegression on dataset:"
    putStrLn . show $ simple
    putStrLn "convergedRegression on dataset:"
    putStrLn . show $ non_robust
    putStrLn "robustFit on dataset:"
    putStrLn . show $ robust
    let filename = "test_robust.png"
    putStrLn $ "Image output is at " ++ filename
    writeFigure PNG filename (800,800) $ testFigure xs ys (simple,non_robust,robust)

randTest :: MonadRandom m => U.Vector Double -> U.Vector Double -> m (EstimatedRelation,EstimatedRelation,EstimatedRelation)
randTest xs ys = do
    robust <- robustFit defaultEstimationParameters xs ys
    let simple = linearRegression xs ys
    let non_robust = converge defaultEstimationParameters xs ys (0.0,0.001) -- simple
    return (simple, non_robust, robust)

test_variances :: IO ()
test_variances = do
    putStrLn "generating random dataset for variance test:"
    let xs = U.fromList [-100..100]
    offsets <- liftM U.fromList $ getNormals 0 10 (U.length xs)
    let ys = U.zipWith (+) xs offsets
    let ab = linearRegression xs ys
    putStrLn $ "estimated fit should be (0,1). It is:" ++ show ab
    let mse = linearRegressionMSE ab xs ys
    putStrLn $ "Calculated MSE of sampled data should be an estimate of 10. it is:" ++ (show . sqrt $ mse)
    let dsts = linearRegressionDistributions ab xs ys
    putStrLn $ "Calculated distributions of the linear fit are linear transformed StudentT distributions with scalings that are estimates of (0.5,1.4777-4). They are:" ++ show dsts

ext: hs

max_stars_repo_name: alpmestan/statistics-linreg

max_issues_repo_name: alpmestan/statistics-linreg

max_forks_repo_name: alpmestan/statistics-linreg


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE Strict #-}
module Pinwheel.Base where

import           Data.Complex

-- {-# INLINE fourierMellin #-}
-- fourierMellin :: (Eq e, RealFloat e) => e -> Int -> Int -> (e, e) -> Complex e
-- fourierMellin sigma angularFreq radialFreq (!x, !y) =
--   if x == 0 && y == 0
--     then 0
--     else (x :+ y) ** (fromIntegral (-angularFreq) :+ 0) *
--          ((x ^ 2 + y ^ 2) :+ 0) **
--          (((fromIntegral angularFreq - 2 + sigma) :+ (-fromIntegral radialFreq)) /
--           2)

{-# INLINE fourierMellin #-}
fourierMellin :: (Eq e, RealFloat e) => e -> Int -> Int -> (e, e) -> Complex e
fourierMellin sigma angularFreq radialFreq (!x, !y) =
  if x == 0 && y == 0
    then 0
    else let r = sqrt $ x ^ 2 + y ^ 2
             theta = atan2 y x
         in ((r :+ 0) ** ((sigma - 1) :+ fromIntegral (-radialFreq))) *
            (cis $ fromIntegral (-angularFreq) * theta)

{-# INLINE fourierMellinInv #-}
fourierMellinInv :: (Eq e, RealFloat e) => e -> Int -> Int -> (e, e) -> Complex e
fourierMellinInv sigma angularFreq radialFreq (!x, !y) =
  if x == 0 && y == 0
    then 0
    else (x :+ y) ** (fromIntegral angularFreq :+ 0) *
         ((x ^ 2 + y ^ 2) :+ 0) **
         (((-(fromIntegral angularFreq) - sigma) :+ fromIntegral radialFreq) /
          2)

{-# INLINE fourierMellinInvPeriod #-}
fourierMellinInvPeriod ::
     (Eq e, RealFloat e) => e -> e -> Int -> Int -> (e, e) -> Complex e
fourierMellinInvPeriod sigma periodEnv angularFreq radialFreq (!x, !y) =
  if x == 0 && y == 0
    then 0
    else (x :+ y) ** (fromIntegral angularFreq :+ 0) *
         ((x ^ 2 + y ^ 2) :+ 0) **
         (((-(fromIntegral angularFreq) - sigma) :+ 2 * pi * fromIntegral radialFreq / (log periodEnv)) /
          2)


text:
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE BangPatterns #-}

{-|
  This module implements PDF instances
  that are missing from random-fu.
  We also include functions to conveniently
  construct values of type SomeDist.
 -}

module Model.PDF where

import Data.Random
import Data.Random.Distribution
import Data.Random.Distribution.Bernoulli hiding ( bernoulli )
import Data.Random.Distribution.Beta hiding ( beta )
import Data.Random.Distribution.Binomial hiding ( binomial )
import Data.Random.Distribution.Categorical hiding ( categorical )
import Data.Random.Distribution.Exponential hiding ( exponential )
import Data.Random.Distribution.Gamma hiding ( gamma )
import Data.Random.Distribution.Normal hiding ( normal )
import Data.Random.Distribution.Poisson hiding ( poisson )
import Data.Random.Distribution.T hiding ( t )
import Data.Random.Distribution.Uniform hiding ( uniform )

import qualified Data.Vector as V

import Model.Types

import Numeric.SpecFunctions

bernoulli :: Double -> SomeDist Bool
bernoulli = SomeDist . Bernoulli

constant :: Eq a => a -> SomeDist a
constant a = categorical [(1, a)]

beta :: Double -> Double -> SomeDist Double
beta a b = SomeDist $ Beta a b

binomial :: Int -> Double -> SomeDist Int
binomial n p = SomeDist $ Binomial n p

normal :: Double -> Double -> SomeDist Double
normal m v = SomeDist $ Normal m v

categorical :: Eq a => [(Double, a)] -> SomeDist a
categorical = SomeDist . fromList

exponential :: Double -> SomeDist Double
exponential = SomeDist . Exp

gamma :: Double -> Double -> SomeDist Double
gamma a b = SomeDist $ Gamma a b

poisson :: Double -> SomeDist Int
poisson = SomeDist . Poisson

uniform :: Double -> Double -> SomeDist Double
uniform l u = SomeDist $ Uniform l u

t :: Integer -> SomeDist Double
t = SomeDist . T

instance PDF (Bernoulli Double) Bool where
  pdf (Bernoulli p) bool = if bool then p else 1 - p
  logPdf d = log . pdf d

instance Eq a => PDF (Categorical Double) a where
  pdf d val = scan l val / (sum $ fst <$> l)
    where scan [] x = 0
          scan ((p, x) : xs) v = if x == v then p else scan xs v
          l = toList d
  logPdf d = log . pdf d

instance PDF Exponential Double where
  pdf (Exp d) x = if x >= 0 then d * exp(-d * x) else 0
  logPdf d = log . pdf d

-- | TODO: Test which parameterization this is.
-- For now, hypothesize that it is the shape-rate parameterization
-- parameterization, where the expectation is (a/b)
-- and the pdf is $ \frac{\beta^\alpha}{\Gamma(\alpha)}  x^{\alpha-1} e^{-\beta x} $.
instance PDF Gamma Double where
  pdf d = exp . logPdf d
  logPdf (Gamma !a !b) !x = a * log b - logGamma a + (a - 1) * log x - b * x

instance PDF (Poisson Double) Int where
  pdf d = exp . logPdf d
  logPdf (Poisson !l) !x = -l + (fromIntegral x * log l) - log (fromIntegral (fact x))

-- | Factorial
fact :: Int -> Int
fact 0 = 1
fact !n = if n < 0 then error "Factorial called on negative" else fact' n 1
  where fact' 1 !m = m
        fact' !m' !m = fact' (m'-1) (m * m')

instance PDF Uniform Double where
  pdf (Uniform !l !u) !x = if x >= l && x <= u then (1 / (u - l)) else 0
  logPdf d = log . pdf d

instance PDF T Double where
  pdf d = exp . logPdf d
  logPdf (T !nu) !x = logGamma ((nud + 1) / 2) - (1/2) * log (nud * pi) - logGamma (nud / 2) + ((-nud-1) / 2) * logGamma (1 + ((x ^ 2) / nud))
    where nud = fromIntegral nu


text:
{- Assignment 1
 - Name: Talha Amjad
 - Date:
 -}
module Assign_1_ExtraCredit where

-- import Data.Complex -- TODO uncomment me to use built-in Complex type
-- see https://www.stackage.org/haddock/lts-8.24/base-4.9.1.0/Data-Complex.html

macid = "amjadt1"


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE ExplicitForAll      #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Elescore.Integration.DB
  ( runDBSource
  ) where

import           ClassyPrelude                          hiding (for, forM, head)
import           Data.DateTime                          (toSqlString)
import           Data.Map                               (elems)
import           Network.HTTP.Client                    (Manager)
import           Pipes
import           Pipes.Concurrent
import qualified Pipes.Prelude                          as P
import           Prelude                                (head)

import           Database.SimpleEventStore
import           Elescore.IdTypes
import           Elescore.Integration.Common.Monitoring
import           Elescore.Integration.Common.Types
import           Elescore.Integration.Common.Utils
import           Elescore.Integration.DB.Client
import           Elescore.Integration.DB.Mapping
import           Elescore.Integration.DB.Monitoring
import           Elescore.Integration.DB.Types
import           Statistics.IQR

runDBSource :: MonadIO m => Store -> ApiKey -> Host -> Manager -> m (Source 'DB)
runDBSource store key host mgr = do
  disP <- disruptions store
  fP <- facilities store
  oEvs <- liftIO (readStream store)

  let oState = replayMkSeed dbObjectMonitor . fmap evPayload $ oEvs
      oIds = keys oState
      oP = objects store oEvs

  (output, input) <- liftIO (spawn unbounded)

  (out1, in1) <- liftIO (spawn unbounded)
  (out2, in2) <- liftIO (spawn unbounded)
  (out3, in3) <- liftIO (spawn unbounded)

  liftIO $ do
    void . forkIO . runEffect $ apiP 10 fetchDisruptedFacilities >-> disP >-> toOutput out1
    void . forkIO . runEffect $ apiP 3600 (fetchFacilities Nothing) >-> fP >-> toOutput (output <> out2)
    void . forkIO . runEffect $ fromInput input >-> P.map evPayload >-> stationFilterP oIds >-> fetchStationP 5 >-> oP >-> toOutput out3

  return Source
     { disruptionEvents = in1
     , facilityEvents = in2
     , objectEvents = in3
     }

  where
    stationFilterP :: MonadIO m => [ObjectId] -> Pipe (FacilityEvent 'DB) ObjectId m ()
    stationFilterP = go
      where
        go s = do
          ev <- await
          case ev of
            FacilityAssignedToObject _ oid -> unless (oid `elem` s) (yield oid >> go (oid : s))
            _                              -> return ()
          go s

    fetchStationP :: Int -> Pipe ObjectId Station IO ()
    fetchStationP delay = forever $ do
      oid <- fromObjectId <$> await
      a <- liftIO $ runReaderT (fetchStation oid) (key, mgr, host)
      either print yield a
      waitSeconds delay

    apiP :: Int -> API a -> Producer a IO ()
    apiP delay action = forever $ do
      a <- liftIO $ runReaderT action (key, mgr, host)
      either print yield a
      waitSeconds delay

disruptions :: MonadIO m => Store -> m (Pipe [Facility] (PersistedEvent (DisruptionEvent 'DB)) IO ())
disruptions store = do
  disEvs <- liftIO (readStream store)

  return $ P.map (fmap mkMDisruption)
    >-> monitorP dbDisruptionMonitor (replayMkSeed dbDisruptionMonitor . fmap evPayload $ disEvs)
    >-> eventStoreP store
    >-> eachBefore (chunkBy (toSqlString . evOccurredOn) disEvs)
    >-> disruptedApiFilter
    >-> P.concat
    >-> filterUnknown

facilities :: MonadIO m => Store -> m (Pipe [Facility] (PersistedEvent (FacilityEvent 'DB)) IO ())
facilities store = do
  fEvs <- liftIO (readStream store)

  return $ P.map (fmap mkMFacility)
    >-> monitorP dbFacilityMonitor (replayMkSeed dbFacilityMonitor . fmap evPayload $ fEvs)
    >-> eventStoreP store
    >-> P.concat
    >-> eachBefore fEvs

objects :: Store -> [PersistedEvent (ObjectEvent 'DB)] -> Pipe Station (PersistedEvent (ObjectEvent 'DB)) IO ()
objects store oEvs =
  P.map mkMObject
    >-> monitorSP dbObjectMonitor (replayMkSeed dbObjectMonitor . fmap evPayload $ oEvs)
    >-> eventStoreP store
    >-> P.concat
    >-> eachBefore oEvs

-- | Filters out FaSta API unknown facility states
filterUnknown :: Pipe (PersistedEvent (DisruptionEvent a)) (PersistedEvent (DisruptionEvent a)) IO ()
filterUnknown = for cat $ \ev ->
  case evPayload ev of
    FacilityDisrupted _ r       -> unless (unknownReason r) (yield ev)
    DisruptionReasonUpdated _ r -> unless (unknownReason r) (yield ev)
    _                           -> yield ev

  where
    unknownReason :: Reason -> Bool
    unknownReason MonitoringNotAvailable = True
    unknownReason MonitoringDisrupted    = True
    unknownReason _                      = False

type Disruptions = MonitorState FacilityId MDisruption

disruptedApiFilter :: MonadIO m => Pipe [PersistedEvent (DisruptionEvent 'DB)] [PersistedEvent (DisruptionEvent 'DB)] m ()
disruptedApiFilter = do
  firstEvents <- await
  let initialState = deriveState mempty firstEvents
      initialSample = singletonSample . fromIntegral . length $ initialState
  go initialState initialSample Nothing

  where
    go :: MonadIO m => Disruptions -> Sample -> Maybe (Disruptions, [PersistedEvent (DisruptionEvent 'DB)]) -> Pipe [PersistedEvent (DisruptionEvent 'DB)] [PersistedEvent (DisruptionEvent 'DB)] m ()
    go previousState sample Nothing = do
      nextEvents <- await
      let nextState = deriveState previousState nextEvents
          numberOfDisruptions = fromIntegral (length nextState)

      if isResidual numberOfDisruptions sample
        then do
        liftIO (putStrLn $ (pack . toSqlString . evOccurredOn $ head nextEvents) <> " Residual detected: " <> tshow (length previousState) <> " -> " <> tshow numberOfDisruptions)
        go previousState sample (Just (nextState, nextEvents))
        else yield nextEvents >> go nextState (addSample numberOfDisruptions sample) Nothing
    go lastKnownGoodState sample (Just (previousState, accumulatedEvents)) = do
      nextEvents <- await
      let nextState = deriveState previousState nextEvents
          numberOfDisruptions = fromIntegral (length nextState)
      if isResidual numberOfDisruptions sample
        then go lastKnownGoodState sample (Just (nextState, accumulatedEvents ++ nextEvents))
        else do
        putStrLn $ (pack . toSqlString . evOccurredOn $ head nextEvents) <> " Back to normal: " <> tshow (length lastKnownGoodState) <> " -> " <> tshow numberOfDisruptions
        yield (compensateApiDisruption lastKnownGoodState nextState (accumulatedEvents ++ nextEvents))
        go nextState (addSample numberOfDisruptions sample) Nothing

    addSample :: Double -> Sample -> Sample
    addSample a = restrictSampleSize 1000 . insertSample' a

    deriveState = foldl' (flip applyEvent)

    applyEvent :: PersistedEvent (DisruptionEvent 'DB) -> Disruptions -> Disruptions
    applyEvent = apply dbDisruptionMonitor . evPayload

    -- | from all happened events, take the latest one that was supposed to be happening
    compensateApiDisruption :: Disruptions -> Disruptions -> [PersistedEvent (DisruptionEvent 'DB)] -> [PersistedEvent (DisruptionEvent 'DB)]
    compensateApiDisruption d1 d2 evs =
      let (_, events) = calculateChanges dbDisruptionMonitor d1 (elems d2)
          reversedEvents = reverse evs
      in foldl' (\acc ev -> acc ++ maybe mempty pure (find ((==) ev . evPayload) reversedEvents)) [] events

chunkBy :: Eq b => (a -> b) -> [a] -> [[a]]
chunkBy _  [] = []
chunkBy f (x:xs) =
  let (_, as, res) = foldl' doChunk (f x, [x], []) xs
  in reverse (reverse as : res)
  where
    doChunk (b, as, ass) a =
      if f a == b
      then (b, a : as, ass)
      else (f a, [a], reverse as : ass)

max_stars_repo_name: Akii/elescore

max_issues_repo_name: Akii/elescore

max_forks_repo_name: Akii/elescore


text:
module Lib where

import Codec.Picture
import System.Directory
import System.IO
import Data.Array.Repa as R
import Data.Array.Repa.Repr.Vector
import Data.Array.Repa.FFTW
import Data.Complex (Complex(..), realPart, magnitude)
import Prelude as P
import Control.Monad

normalize :: Monad m => Array U DIM2 Double -> m (Array D DIM2 Double)
normalize arr = do
    minn <- foldAllP min 1 arr
    let arrmin = R.map (+ (-minn)) arr
    maxx <- foldAllP max 0 arrmin
    let arrnorm = R.map (/maxx) arrmin
    return arrnorm

myfilter tresh v
    | v < tresh = fromIntegral v * 0.6
    | otherwise = 0

rowFilter :: Monad m => Array D DIM1 Double -> m (Array D DIM1 Double)
rowFilter row = do
    rowComplex <- computeP $ R.map (\x -> x :+ 0 ) row
    let fftF = fft $ rowComplex
        (Z :. w) = extent fftF
        barier = round (fromIntegral w/2)
        filtered = fromListUnboxed (Z :. w) $ (P.map (myfilter barier) [1..w])
        fftFilt = R.zipWith (*) fftF filtered
    ifftF <- fmap ifft $ computeP fftFilt
    return $ R.map realPart ifftF

getRow :: Array D DIM2 Double -> Int -> Array D DIM1 Double
getRow array n = slice array (Any :. n :. All)

mapRows :: Monad m => (Array D DIM1 Double -> m (Array D DIM1 Double)) -> Array D DIM2 Double -> m (Array U DIM2 Double)
mapRows func array = do
    let (Z :. w :. h) = extent array
    rows <- mapM (\num -> do
        let row = getRow array num
        func row) [0,1..(w-1)]
    let hugeRow = toList $ foldr1 append rows
    return $ fromListUnboxed (Z :. w :. h) hugeRow

reconstruct :: Array D DIM2 Double -> Double -> Int -> IO ()
reconstruct img p orgW = do
    let (Z :. w :. h) = extent img
        angleStep = pi / fromIntegral h
        anglesList = takeWhile (<pi) [a * angleStep | a <- [0..]]
        orgWNum = fromIntegral orgW

        listOfIndicies (x, y) = let
            list = P.map (\a -> round $ (x * sin a + y * cos a)/p + orgWNum/2) anglesList
            list_zip = zip list [0..(h-1)]
            in [(a,b) | (a,b) <- list_zip, a >= 0, a < w]

        render (x, y) = let
            list = listOfIndicies (x, y)
            pixelList = P.map (\(p, h) -> img ! (Z :. p :. h)) list
            pixelSum = sum pixelList
            avg = pixelSum / (fromIntegral $ length pixelList)
            in if avg > 0 then avg else 0

    let imageIndicies = [orgWNum/(-2)..orgWNum/2-1]
        img' = fromListUnboxed (Z :. orgW :. orgW) (P.map render [(a,b) | a <- imageIndicies , b <- imageIndicies])
    img <- normalize img'
    writePng "res/reconstruct.png" $ generateImage (\x y -> dToPx (img ! (Z :. x :. y))) orgW orgW

getY a p r x = (x, round ((p - (fromIntegral (x - r) * (cos a))) / (sin a)) + r)
getX a p r y = (round ((p - (fromIntegral (y - r) * (sin a))) / (cos a)) + r, y)

filterCoords r (x, y) = ((fromIntegral x - r)**2 + (fromIntegral y - r)**2) < r**2

getLineAvg :: Double -> Double -> Int -> (Double -> Double -> Int -> Int -> (Int, Int)) -> Array U DIM2 Double -> Double
getLineAvg a p w getCoord img =
    let pixelList = P.map (getCoord a p (div w 2)) [0..w-1]
        r = (fromIntegral w) / 2
        pixelList' = filter (filterCoords r) pixelList
        len = length pixelList'
        ret = (sum $ P.map (\(x, y) -> img ! (Z :. x :. y)) pixelList') / fromIntegral len
    in if isNaN ret then 0 else ret

getDetectorValue :: Double -> Double -> Int -> Array U DIM2 Double -> Double
getDetectorValue a p w img
    | a < pi/4 || a > (3 * pi)/4 = getLineAvg a p w getX img
    | otherwise = getLineAvg a p w getY img

processImage :: String -> Int -> Int -> Int -> IO ()
processImage fname nsteps nrays opening' = do
    (w, h, img) <- imageToArray fname
    putStrLn "Image loaded"
    let opening = (fromIntegral opening'/ 180) * pi
        r = fromIntegral w/2
        tresh = sin (opening/2) * r
        p = (sin (opening/2) * r) * 2 / (fromIntegral nrays)
        step = pi / fromIntegral nsteps
        angle = takeWhile (<pi) [x * step | x <- [0..]]
    putStrLn "Calculating projections"
    let listOfP = P.map (\v -> ((fromIntegral v) * p) - tresh) [0..nrays-1]
        rest = P.map (\a -> P.map (\p -> getDetectorValue a p w img) listOfP) angle
        trans = fromListUnboxed (Z :. length rest :. length (head rest)) $ concat $ rest
        trans' = transpose trans
    x <- computeUnboxedP trans'
    nofilter <- normalize x
    arraySaveToImage nofilter "nofilter"
    result' <- mapRows rowFilter nofilter
    result <- normalize result'
    arraySaveToImage result "filter"
    reconstruct result p w

dToPx x
    | x < 0 = PixelYA8 0 255
    | x > 1 = PixelYA8 255 255
    | otherwise = PixelYA8 (round (x * 255)) 255

arraySaveToImage :: Array D DIM2 Double -> String -> IO ()
arraySaveToImage arr fname = do
    let (Z :. w :. h) = extent arr
    writePng ("res/" P.++ fname P.++ ".png") $ generateImage (\x y -> dToPx $ arr ! (Z :. x :. y)) w h

rgbToGrey :: PixelRGB8 -> Double
rgbToGrey (PixelRGB8 r g b) =
  let rgb = P.zipWith (*) [0.2126, 0.7152, 0.0722] $ P.map (fromIntegral) [r, g, b]
  in fromInteger . round . sum $ rgb :: Double

imageToArray :: String -> IO (Int, Int, Array U DIM2 Double)
imageToArray fname = do
    Right imgRGB <- readImage fname
    let (imgGray@(Image w h _)) = convertRGB8 imgRGB
        arr = fromFunction (Z :. w :. h) (\(Z :. x :. y) -> rgbToGrey $ pixelAt imgGray x y)
    img <- computeUnboxedP arr
    return (w, h, img)


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE FlexibleContexts #-}

module LinAlg where

import qualified Numeric.LinearAlgebra as L

import Data.Complex
import Data.List

class Vector v where
    -- vector addition
    (<+>) :: v -> v -> v
    -- scalar multiplication
    (<**>) :: Complex Double -> v -> v
    -- additive inverse
    vneg :: v -> v
    (<~>) :: v -> v -> v
    -- vector subtraction for convenience
    u <~> v = u <+> vneg v

class Vector v => Hilbert v where
    -- Inner product
    (<.>) :: v -> v -> Complex Double
    -- Vector norm, default definition generates norm
    -- from inner product
    norm :: v -> Double
    norm v = sqrt . magnitude $ v <.> v

-- Data structure representing a generic linear operator
-- in some representation
newtype Operator r = Operator {
    repr :: r
}

instance Functor Operator where
    fmap f (Operator r) = Operator $ f r

instance Applicative Operator where
    pure = Operator
    (Operator f) <*> op = f <$> op

-- Type alias for matrix representation of operator
type MatOp = Operator (L.Matrix (Complex Double))

-- Allow operators to be composed
instance Semigroup r => Semigroup (Operator r) where
    (Operator r1) <> (Operator r2) = Operator (r1 <> r2)

-- Operators are vectors
instance Vector MatOp where
    o1 <+> o2 = Operator (repr o1 + repr o2)
    a <**> o = L.cmap (*a) <$> o
    vneg o = negate <$> o

-- Trace of an operator
trace :: MatOp -> Complex Double
trace = L.sumElements . L.takeDiag . repr

max_stars_repo_name: Jvanrhijn/optiqs

max_issues_repo_name: Jvanrhijn/optiqs

max_forks_repo_name: Jvanrhijn/optiqs


text:
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE GADTs #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}

module UntypedNet where

import Data.Functor
import Data.Maybe
import Data.List (foldl')
import Numeric.LinearAlgebra (Matrix, Vector, randomVector, uniformSample, RandDist(Uniform), (#>), tr, outer, scale, norm_2, vector)
import Control.Monad.Random (MonadRandom, getRandom, replicateM)

-- Wx + b
data Weights = W { wBiases :: !(Vector Double)  -- n
                 , wNodes :: !(Matrix Double)  -- n x m
                 }  -- "m to n" layer

-- define a network data type
-- can build a network like : ih :&~ hh :&~ O ho
data Network :: * where
    O :: !Weights -> Network
    (:&~) :: !Weights -> !Network -> Network

infixr 5 :&~

{- equivalent to:
data Network = O Weights
             | Weights :&~ Network
-}

randomWeights :: MonadRandom m => Int -> Int -> m Weights
randomWeights i o = do
  seed1 :: Int <- getRandom  -- generate random seed
  seed2 :: Int <- getRandom
      -- randomVector :: Seed -> RandDist -> Int -> VectorDouble
      -- obtains a vector of random elements of provided size
  let wB = randomVector seed1 Uniform o * 2 - 1
      -- uniformSample :: Seed -> Int -> [(Double, Double)] -> Matrix Double
      wN = uniformSample seed2 o (replicate i (-1, 1))  -- o x i -sized matrix
  return $ W wB wN

-- build a random network composed of layers
-- input layer size -> [hidden layer sizes] -> output layer size
randomNet :: MonadRandom m => Int -> [Int] -> Int -> m Network
randomNet i [] o = O <$> randomWeights i o  -- final fully-connected layer
randomNet i (h:hs) o = (:&~) <$> randomWeights i h <*> randomNet h hs o

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

-- derivative function of logistic
logistic' :: Floating a => a -> a
logistic' x = logix * (1 - logix)
  where
    logix = logistic x

-- matrix-vector multiplication
-- (#>) :: Numeric t => Matrix t -> Vector t -> Vector t
runLayer :: Weights -> Vector Double -> Vector Double
runLayer (W wB wN) v = wB + wN #> v

-- feedforward the network with every activaciton function as sigmoid
runNet :: Network -> Vector Double -> Vector Double
runNet (O w) !v = logistic (runLayer w v)
runNet (w :&~ n') !v = let v' = logistic (runLayer w v)
                       in runNet n' v'

train :: Double  -- learning rate
      -> Vector Double  -- input vector
      -> Vector Double  -- target vector
      -> Network  -- network to train
      -> Network
train rate x0 target = fst . go x0
  where
    go :: Vector Double  -- input vector
       -> Network  -- network to train
       -> (Network, Vector Double)  -- (new network layer, back-propping derivative chain)
    -- handle the output layer
    go !x (O w@(W wB wN))
      = let y = runLayer w x
            o = logistic y
            -- the gradient (dE / dy)
            dEdy = logistic' y * (o - target)
            -- bias and node update
            wB' = wB - scale rate dEdy  -- multiply by scalar 'rate'
            wN' = wN - scale rate (dEdy `outer` x)
            w' = W wB' wN'
            -- bundle of derivatives for next step
            dWs = tr wN #> dEdy
         in (O w', dWs)
    -- handle the inner layers
    go !x (w@(W wB wN) :&~ n)
        = let y = runLayer w x
              o = logistic y
              -- get dWs', bundle of derivatives from rest of the net
              (n', dWs') = go o n
              -- the gradient (how much y affects the error)
              dEdy = logistic' y * dWs'
              -- new bias weights and node weights
              wB' = wB - scale rate dEdy
              wN' = wN - scale rate (dEdy `outer` x)
              w' = W wB' wN'
              -- bundle of derivatives for next step
              dWs = tr wN #> dEdy
          in (w' :&~ n', dWs)

netTest :: MonadRandom m => Double -> Int -> m String
netTest rate n = do
    inps <- replicateM n $ do
      s <- getRandom
      return $ randomVector s Uniform 2 * 2 - 1
    let outs = flip map inps $ \v ->
                 if v `inCircle` (fromRational 0.33, 0.33)
                      || v `inCircle` (fromRational (-0.33), 0.33)
                   then fromRational 1
                   else fromRational 0
    net0 <- randomNet 2 [16,8] 1
    let trained = foldl' trainEach net0 (zip inps outs)
          where
            trainEach :: Network -> (Vector Double, Vector Double) -> Network
            trainEach nt (i, o) = train rate i o nt

        outMat = [ [ render (norm_2 (runNet trained (vector [x / 25 - 1,y / 10 - 1])))
                   | x <- [0..50] ]
                 | y <- [0..20] ]
        render r | r <= 0.2  = ' '
                 | r <= 0.4  = '.'
                 | r <= 0.6  = '-'
                 | r <= 0.8  = '='
                 | otherwise = '#'

    return $ unlines outMat
  where
    inCircle :: Vector Double -> (Vector Double, Double) -> Bool
    v `inCircle` (o, r) = norm_2 (v - o) <= r

(!!?) :: [a] -> Int -> Maybe a
xs !!? i = listToMaybe (drop i xs)

max_stars_repo_name: deNsuh/hasml

max_issues_repo_name: deNsuh/hasml

max_forks_repo_name: deNsuh/hasml


text:
{-# LANGUAGE TupleSections #-}

import Control.Monad.IO.Class (liftIO)
import Data.Complex (Complex ((:+)))
import Data.Functor.Classes (liftEq2)
import Foreign (Storable)
import Streamly.Internal.Data.Stream.IsStream (SerialT)
import Test.Hspec.Core.Spec (SpecM)
import Test.Hspec.QuickCheck (prop)
import Test.QuickCheck
    (elements, chooseInt, choose, forAll, Property, vectorOf)
import Test.QuickCheck.Monadic (monadicIO, assert)

import qualified Data.Map.Strict as Map
import qualified Data.Set as Set
import qualified Data.Vector as V
import qualified Statistics.Sample.Powers as STAT
import qualified Statistics.Transform as STAT
import qualified Streamly.Internal.Data.Array.Foreign.Mut as MA
import qualified Streamly.Internal.Data.Array.Foreign.Type as Array
import qualified Streamly.Internal.Data.Fold as Fold
import qualified Streamly.Internal.Data.Ring.Foreign as Ring
import qualified Streamly.Internal.Data.Stream.IsStream as Stream
import qualified Streamly.Prelude as S

import Prelude hiding (sum, maximum, minimum)

import Streamly.Statistics
import Test.Hspec

tolerance :: Double
tolerance = 0.00001

validate :: Double -> Bool
validate delta  = delta < tolerance

jackKnifeInput :: [Double]
jackKnifeInput = [1.0::Double, 2.0, 3.0, 4.0]

jackMeanRes :: [Double]
jackMeanRes = [3.0, 2.6666666666666665, 2.3333333333333335, 2.0]

jackVarianceRes :: [Double]
jackVarianceRes =
    [ 0.6666666666666661
    , 1.5555555555555554
    , 1.5555555555555545
    , 0.666666666666667
    ]

jackStdDevRes :: [Double]
jackStdDevRes =
    [ 0.8164965809277257
    , 1.247219128924647
    , 1.2472191289246466
    , 0.8164965809277263
    ]

testDistributions
    :: (STAT.Powers -> Double)
    -> Fold.Fold IO (Double, Maybe Double) Double
    -> Property
testDistributions func fld =
    forAll (chooseInt (1, 1000)) $ \list_length ->
        forAll (vectorOf list_length (choose (-50.0 :: Double, 100.0)))
            $ \ls ->
                monadicIO $ do
                let var2 = func . STAT.powers 2 $ V.fromList ls
                    strm = S.fromList ls
                var1 <-
                    liftIO $ S.fold (Ring.slidingWindow list_length fld) strm
                assert (validate $ abs (var1 - var2))

testVariance :: Property
testVariance = testDistributions STAT.variance variance

testStdDev :: Property
testStdDev = testDistributions STAT.stdDev stdDev

testFuncMD ::
    Fold.Fold IO ((Double, Maybe Double), IO (MA.Array Double)) Double -> Spec
testFuncMD f = do
                let c = S.fromList [10.0, 11.0, 12.0, 14.0]
                a1 <- runIO $ S.fold (Ring.slidingWindowWith 2 f) c
                a2 <- runIO $ S.fold (Ring.slidingWindowWith 3 f) c
                a3 <- runIO $ S.fold (Ring.slidingWindowWith 4 f) c
                it ("MD should be 1.0 , 1.1111111111111114 , 1.25 but actual is "
                    ++ show a1 ++ " " ++ show a2 ++ " " ++ show a3)
                    (  validate (abs (a1 - 1.0))
                    && validate (abs (a2 - 1.1111111))
                    && validate (abs (a3 - 1.25))
                    )

testFuncKurt :: Spec
testFuncKurt = do
    let c = S.fromList
            [21.3 :: Double, 38.4, 12.7, 41.6]
    krt <- runIO $ S.fold (Ring.slidingWindow 4 kurtosis) c
    it ( "kurtosis should be 1.2762447351370185 Actual is " ++
        show krt
        )

        (validate $ abs (krt - 1.2762447))

testJackKnife :: (Show a, Eq a, Storable a) =>
       (Array.Array a -> SerialT (SpecM ()) a)
    -> [a]
    -> [a]
    -> Spec
testJackKnife f ls expRes = do
    let arr = Array.fromList ls
    res <- Stream.toList $ f arr
    it ("testJackKnife result should be ="
        ++ show expRes
        ++ " Actual is = " ++show res
        )
        (res == expRes)

testFuncHistogram :: Spec
testFuncHistogram = do
    let strm = S.fromList [1..15]
    res <- runIO $
        S.fold (histogram (binOffsetSize (0::Int) (3::Int))) strm
    let expected = Map.fromList
                    [ (0::Int, 2::Int)
                    , (1, 3)
                    , (2, 3)
                    , (3, 3)
                    , (4, 3)
                    , (5, 1)
                    ]

    it ("Map should be = "
        ++ show expected
        ++ " Actual is = "
        ++ show res) (expected == res)

testFuncbinFromSizeN :: Int -> Int -> Int -> Int -> HistBin Int -> SpecWith (Arg Bool)
testFuncbinFromSizeN low binSize nbins x exp0 = do
    let res = binFromSizeN low binSize nbins x
    it ("Bin should be = "
        ++ show exp0
        ++ " Actual is = "
        ++ show res) (res == exp0)

testFuncbinFromToN :: Int -> Int -> Int -> Int -> HistBin Int -> SpecWith ()
testFuncbinFromToN low high n x exp0 = do
    let res = binFromToN low high n x
    it ("Bin should be = "
        ++ show exp0
        ++ " Actual is = "
        ++ show res) (res == exp0)

testFrequency :: [Int] -> Map.Map Int Int -> Spec
testFrequency inputList result = do
    freq <- S.fold frequency $ S.fromList inputList
    it ("Frequency " ++ show freq) $ liftEq2 (==) (==) freq result

testMode :: [Int] -> Maybe (Int, Int) -> Spec
testMode inputList res = do
    mode0 <- S.fold mode $ S.fromList inputList
    it ("Mode " ++ show mode0) $ mode0 == res

testFFT :: Property
testFFT = do
    let lengths = [2, 4, 8, 16]
    forAll (elements lengths) $ \list_length ->
        forAll (vectorOf list_length (choose (-50.0 :: Double, 100.0)))
            $ \ls ->
                monadicIO $ do
                    let tc = map (\x -> x :+ 0) ls
                    let vr = V.toList (STAT.fft (V.fromList tc)
                                        :: V.Vector STAT.CD)
                    marr <- MA.fromList tc
                    fft marr
                    res <- MA.toList marr
                    assert (vr == res)

sampleList :: [Double]
sampleList = [1.0, 2.0, 3.0, 4.0, 5.0]

testResample :: [Double] -> Spec
testResample sample = do
    let sampleArr = Array.fromList sample
        sampleSet = Set.fromList sample
    resampleList <- runIO $ S.toList $ S.unfold resample sampleArr
    let resampleSet = Set.fromList resampleList
        sub = Set.isSubsetOf resampleSet sampleSet
    -- XXX We should not use dynamic output in test description
    it ("resample " ++ show resampleList)
       (Prelude.length resampleList == Array.length sampleArr && sub)

testFoldResamples :: Int -> [Double] -> Spec
testFoldResamples n sample = do
    let arr = Array.fromList sample
    a <- runIO $ S.toList $ foldResamples n arr Fold.mean
    -- XXX We should not use dynamic output in test description
    it ("foldResamples " ++ show a) (Prelude.length a == n)

main :: IO ()
main = hspec $ do
    describe "Numerical stability while streaming" $ do
        let numElem = 80000
            winSize = 800
            testCaseChunk = [9007199254740992, 1, 1.0 :: Double,
                                9007199254740992, 1, 1, 1, 9007199254740992]
            testCase = take numElem $ cycle testCaseChunk
            deviationLimit = 1
            testFunc f = do
                let c = S.fromList testCase
                a <- runIO $ S.fold (Ring.slidingWindow winSize f) c
                b <- runIO $ S.fold f $ S.drop (numElem - winSize)
                        $ S.map (, Nothing) c
                let c1 = a - b
                it ("should not deviate more than " ++ show deviationLimit)
                    $ c1 >= -1 * deviationLimit && c1 <= deviationLimit

        describe "Sum" $ testFunc sum
        describe "mean" $ testFunc mean
        describe "welfordMean" $ testFunc welfordMean

    describe "Correctness" $ do
        let winSize = 3
            testCase1 = [31, 41, 59, 26, 53, 58, 97] :: [Double]
            testCase2 = replicate 5 1.0 ++ [7.0]

            testFunc tc f sI sW = do
                let c = S.fromList tc
                a <- runIO $ S.toList $ S.postscan f $ S.map (, Nothing) c
                b <- runIO $ S.toList $ S.postscan
                        (Ring.slidingWindow winSize f) c
                it "Infinite" $ a  == sI
                it ("Finite " ++ show winSize) $ b == sW

        -- Resampling
        describe "JackKnife Mean" $
            testJackKnife jackKnifeMean jackKnifeInput jackMeanRes
        describe "JackKnife Variance" $ do
            testJackKnife jackKnifeVariance jackKnifeInput jackVarianceRes
        describe "JackKnife StdDev" $
            testJackKnife jackKnifeStdDev jackKnifeInput jackStdDevRes

        describe "resample" $ do
            testResample sampleList
        describe "foldResamples 4" $ do
            testFoldResamples 4 sampleList
        describe "foldResamples 6" $ do
            testFoldResamples 6 sampleList

        -- Spread/Mean
        describe "MD" $ testFuncMD md
        describe "Kurt" testFuncKurt
        prop "fft" testFFT
        describe "minimum" $ do
            let scanInf = [31, 31, 31, 26, 26, 26, 26] :: [Double]
                scanWin = [31, 31, 31, 26, 26, 26, 53] :: [Double]
            testFunc testCase1 minimum scanInf scanWin
        describe "maximum" $ do
            let scanInf = [31, 41, 59, 59, 59, 59, 97] :: [Double]
                scanWin = [31, 41, 59, 59, 59, 58, 97] :: [Double]
            testFunc testCase1 maximum scanInf scanWin
        describe "range" $ do
            let scanInf = [0, 10, 28, 33, 33, 33, 71] :: [Double]
                scanWin = [0, 10, 28, 33, 33, 32, 44] :: [Double]
            testFunc testCase1 range scanInf scanWin
        describe "sum" $ do
            let scanInf = [1, 2, 3, 4, 5, 12] :: [Double]
                scanWin = [1, 2, 3, 3, 3, 9] :: [Double]
            testFunc testCase2 sum scanInf scanWin
        describe "mean" $ do
            let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]
                scanWin = [1, 1, 1, 1, 1, 3] :: [Double]
            testFunc testCase2 mean scanInf scanWin
        describe "welfordMean" $ do
            let scanInf = [1, 1, 1, 1, 1, 2] :: [Double]
                scanWin = [1, 1, 1, 1, 1, 3] :: [Double]
            testFunc testCase2 welfordMean scanInf scanWin

        -- Probability Distribution
        describe "frequency"
            $ testFrequency
                [1::Int, 1, 2, 3, 3, 3]
                (Map.fromList [(1, 2), (2, 1), (3, 3)])
        describe "Mode" $ testMode [1::Int, 1, 2, 3, 3, 3] (Just (3, 3))
        describe "Mode Empty " $ testMode ([]::[Int]) Nothing

        describe "histogram" testFuncHistogram
        describe "binFromSizeN AboveRange" $
            testFuncbinFromSizeN (0::Int) 2 10 55 AboveRange
        describe "binFromSizeN BelowRange" $
            testFuncbinFromSizeN (0::Int) 2 10 (-1) BelowRange
        describe "binFromSizeN InRange" $
            testFuncbinFromSizeN (0::Int) 2 10 19 (InRange 9)
        describe "binFromSizeN AboveRange" $
            testFuncbinFromSizeN (0::Int) 2 10 20 AboveRange
        describe "binFromToN AboveRange" $
            testFuncbinFromToN (0::Int) 49 10 55 AboveRange
        describe "binFromToN BelowRange" $
            testFuncbinFromToN (0::Int) 49 10 (-1) BelowRange
        describe "binFromToN InRange"    $
            testFuncbinFromToN (0::Int) 49 10 19 (InRange 3)
        describe "binFromToN AboveRange" $
            testFuncbinFromToN (0::Int) 50 10 20 (InRange 4)
        prop "variance" testVariance
        prop "stdDev" testStdDev

max_stars_repo_name: composewell/streamly-statistics

max_issues_repo_name: composewell/streamly-statistics

max_forks_repo_name: composewell/streamly-statistics


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE BlockArguments #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE BangPatterns #-}
-- | A music theory library for just intonation and other mathematically pure ideas.
module Boopadoop
  (module Boopadoop
  ,module Boopadoop.Diagram
  ,module Boopadoop.Rhythm
  ,module Boopadoop.Interval
  ,module Boopadoop.Discrete
  ) where

import qualified Data.WAVE as WAVE
import Control.Monad.ST
import Control.Monad
import Control.Applicative
import Boopadoop.Diagram
import Boopadoop.Rhythm
import Boopadoop.Interval
import Boopadoop.Discrete
import Data.List
import Data.Int
import Data.Complex
import qualified Data.Vector.Fixed as FV
import qualified Data.Primitive.ByteArray as BA
import qualified Data.Vector.Unboxed as Vector

-- | A 'Waveform' is a function (of time) that we can later sample.
newtype Waveform t a = Waveform
  {sample :: t -> a -- ^ 'sample' the 'Waveform' at a specified time
  }

instance Functor (Waveform t) where
  fmap f w = sampleFrom $ f . sample w

instance SummaryChar (Waveform t a) where
  sumUp _ = '~'

-- | A 'Double' valued wave with time also in terms of 'Double'.
-- This models a real-valued waveform which typically has values in @[-1,1]@ and
-- is typically supported on either the entire real line ('sinWave') or on a compact subset ('compactWave')
type DWave = Waveform Double Double

-- | Show a waveform by pretty printing some of the actual waveform in dot matrix form.
instance Show (Waveform Double Double) where
  show w = intercalate "\n" . transpose $ map sampleToString waveSamples
    where
      sampleToString k = if k <= quantLevel && k >= -quantLevel
        then replicate (quantLevel - k) '.' ++ "x" ++ replicate (quantLevel + k) '.'
        else let m = "k = " ++ show k in m ++ replicate (quantLevel * 2 + 1 - length m) ' '
      waveSamples = map (floor . (* realToFrac quantLevel) . sample w . (/sampleRate)) [0 .. 115]
      quantLevel = 15 :: Int
      sampleRate = 6400

instance Show (Waveform Double Discrete) where
  show w = intercalate "\n" . transpose $ map sampleToString waveSamples
    where
      sampleToString k = if k <= quantLevel && k >= -quantLevel
        then replicate (quantLevel - k) '.' ++ "x" ++ replicate (quantLevel + k) '.'
        else let m = "k = " ++ show k in m ++ replicate (quantLevel * 2 + 1 - length m) ' '
      waveSamples = map ((`quotRoundUp` (1 + (discFactor `quot` quantLevel))) . fromIntegral . unDiscrete . sample w . (/sampleRate)) [0 .. 115]
      quantLevel = 15 :: Int
      sampleRate = 6400

instance Show (Waveform Tick Discrete) where
  show w = intercalate "\n" . transpose $ map sampleToString waveSamples
    where
      sampleToString k = if k <= quantLevel && k >= -quantLevel
        then replicate (quantLevel - k) '.' ++ "x" ++ replicate (quantLevel + k) '.'
        else let m = "k = " ++ show k in m ++ replicate (quantLevel * 2 + 1 - length m) ' '
      waveSamples = map ((`quotRoundUp` (1 + (discFactor `quot` quantLevel))) . fromIntegral . unDiscrete . sample (skipTicks 1 w)) [0 .. 115]
      quantLevel = 15 :: Int

instance Show (Waveform Tick Double) where
  show w = intercalate "\n" . transpose $ map sampleToString waveSamples
    where
      sampleToString k = if k <= quantLevel && k >= -quantLevel
        then replicate (quantLevel - k) '.' ++ "x" ++ replicate (quantLevel + k) '.'
        else let m = "k = " ++ show k in m ++ replicate (quantLevel * 2 + 1 - length m) ' '
      waveSamples = map (floor . (* realToFrac quantLevel) . sample (skipTicks 1 w)) [0 .. 115]
      quantLevel = 15 :: Int

-- | A version of @'quot'@ that rounds away from zero instead of towards it.
quotRoundUp :: Int -> Int -> Int
quotRoundUp a b = if a `mod` b == 0 then a `quot` b else (signum a * signum b) + (a `quot` b)

-- | Build a 'Waveform' by sampling the given function.
sampleFrom :: (t -> a) -> Waveform t a
sampleFrom f = Waveform $ \t -> {-t `seq`-} f t

-- | Sample a 'Waveform' at specified time. @'sampleAt' = 'flip' 'sample'@
sampleAt :: t -> Waveform t a -> a
sampleAt = flip sample

-- | Pure sine wave of the given frequency
--sinWave :: Floating a => a -> Waveform a a
sinWave :: Double -> DWave
sinWave f = if sinHack -- speedy fast sin hack from https://www.gamedev.net/forums/topic/621589-extremely-fast-sin-approximation/ due to nightcracker
  then sampleFrom $ \t -> let
    x = pi * ((2 * f * t) - fromIntegral k)
    x2 = x * x
    k = spookyFastTruncate (2 * f * t)
    in (if odd k then negate else id) (x * (c + x2 * (b + a * x2)))
  else sampleFrom $ \t -> sin (freq * t)
  where
    sinHack = True
    !freq = 2 * pi * f
    a = 0.00735246819687011731341356165096815
    b = -0.16528911397014738207016302002888890
    c = 0.99969198629596757779830113868360584

-- | Pure sine wave of the given frequency and initial phase
sinWaveWithPhase :: Floating a => a -> a -> Waveform a a
sinWaveWithPhase f phi0 = sampleFrom $ \t -> let !freq = 2 * pi * f in sin (freq * t + phi0)

-- | Sine wave that is optimized to store only a small @'CompactWavetable'@.
-- The period is given in ticks because non-tick-multiple period sin waves will have aliasing under this optimization!
fastSin :: Tick -> Wavetable
fastSin per = exploitPeriodicity per $ tickTable (fromIntegral per) $ discretize $ sinWave 1

-- | @'compactWave' (l,h)@ is a wave which is @'True'@ on @[l,h)@ and @'False'@ elsewhere
compactWave :: (Ord t,Num t) => (t,t) -> Waveform t Bool
compactWave (low,high) = sampleFrom $ \t -> t >= low && t < high

-- | @'muting' 'True'@ is @'id'@ while @'muting' 'False'@ is @'const' 0@.
muting :: Num a => Bool -> a -> a
muting b s = if b then s else 0

-- | Modulate one wave with another according to the given function pointwise.
-- This means you can't implement 'phaseModulate' using only this combinator because phase modulation
-- requires information about the target wave at times other than the current time.
modulate :: (a -> b -> c) -> Waveform t a -> Waveform t b -> Waveform t c
modulate f a b = sampleFrom $ \t -> f (sample a t) (sample b t)

-- | Modulate the amplitude of one wave with another. This is simply pointwise multiplication:
-- @
--  'amplitudeModulate' = 'modulate' ('*')
-- @
amplitudeModulate :: Num a => Waveform t a -> Waveform t a -> Waveform t a
amplitudeModulate = modulate (*)

-- | Control rate of speed of one waveform with another
freqModulate :: Waveform t' p -> Waveform p a -> Waveform t' a
freqModulate phaseSelect w = sampleFrom $ \t -> sample w (sample phaseSelect t)

-- | First argument is valued in sample time steps, so @1/stdtr@ is a good starting point.
emuVCO :: Waveform Tick Double -> Waveform Double a -> [a]
emuVCO freqSelect w = go 0 0
  where
    go !i sampPoint = let skip = sample freqSelect i in sample w sampPoint : go (i + 1) (sampPoint + skip)

emuVCO' :: [Double] -> Waveform Double a -> [a]
emuVCO' fs w = go fs 0
  where
    go freqSelect !sampPoint = let (skip:ss) = freqSelect in sample w sampPoint : go ss (sampPoint + skip)

niceVCO :: Double -> [Double] -> Waveform Double a -> [a]
niceVCO amp fs = emuVCO' (fmap (\f -> (amp ** f)/stdtr) fs)

oscAboveZero :: Double -> Double -> Double
oscAboveZero amp = oscAbout (amp/2) amp

oscAbout :: Double -> Double -> Double -> Double
oscAbout z0 amp z = (amp * z) + z0

-- | Piecewise linear reconstruction (interpolation) of a sampled signal
linearResample :: Double -> Waveform Tick Double -> Waveform Double Double
linearResample sampleRate w = sampleFrom $ \t -> do
  let t' = t * sampleRate
  let (low,high) = (floor t',ceiling t')
  (fromIntegral high - t') * sample w low + (t' - fromIntegral low) * sample w high

recordStream :: [a] -> Waveform Tick a
recordStream xs = sampleFrom ((xs!!) . fromIntegral)

-- | The output wavetables here have @sample (toConvolutionTableStream ds !! t) k = ds !! (t - k)@ for @k < windowSize@
toConvolutionTableStream :: Tick -> [Discrete] -> [Wavetable]
toConvolutionTableStream windowSize xs = fmap toWavetable . scanl stepSt (startTable xs) $ xs
  where
    toWavetable :: (BA.ByteArray,Int) -> Wavetable
    toWavetable (v,b) = sampleFrom $ \k -> Discrete $ BA.indexByteArray v ((b - k) `mod` windowSize)
    startTable :: [Discrete] -> (BA.ByteArray,Int)
    startTable ds = (BA.byteArrayFromListN windowSize (fmap unDiscrete $ take windowSize ds),0)
    -- Table, breakpoint
    stepSt :: (BA.ByteArray, Int) -> Discrete -> (BA.ByteArray, Int)
    stepSt (!t,!b) x = runST $ do
      v <- BA.unsafeThawByteArray t
      BA.writeByteArray v (b `mod` windowSize) (unDiscrete x)
      v' <- BA.unsafeFreezeByteArray v
      pure (v',(b+1) `mod` windowSize)

streamWavetable :: Waveform Tick a -> [a]
streamWavetable w = map (sample w) [0..]

memorizeStream :: Int -> [Discrete] -> Wavetable
memorizeStream len = fromCompact . CompactWavetable . Vector.fromList . map unDiscrete . take len

-- | Very slow
integrateWave :: Num a => Waveform Tick a -> Waveform Tick a
integrateWave w = sampleFrom $ \t -> sum . map (sample w) $ [0 .. t]

-- | Modulate the phase of one wave with another. Used in synthesis.
-- @
--  'phaseModulate' beta ('setVolume' 0.2 $ 'sinWave' 'concertA') ('setVolume' 0.38 $ 'triWave' 'concertA')
-- @
-- (try beta=0.0005)
phaseModulate :: Num t
              => t -- ^ Tuning parameter. Modulation signal is @'amplitudeModulate'@d by @('const' beta)@
              -> Waveform t t -- ^ Modulation signal. Outputs the phase shift to apply
              -> Waveform t a -- ^ Target wave to be modulated
              -> Waveform t a
phaseModulate beta modulation target = sampleFrom $ \t -> sample target (t + beta * sample modulation t)

-- | Smoothly transition to playing a wave back at a different speed after some time
changeSpeed :: (Ord a,Fractional a) => a -> a -> a -> Waveform a a -> Waveform a a
changeSpeed startTime lerpTime newSpeed wave = sampleFrom $ \t -> sample wave $ if t < startTime
  then t
  else if t > startTime + lerpTime
    then startTime + newSpeed * t
    -- Lerp between sampling at 1 and sampling at newSpeed
    else startTime + (1 + ((t - startTime)/lerpTime) * (newSpeed - 1)) * t

-- | Play several waves on top of each other, normalizing so that e.g. playing three notes together doesn't triple the volume.
balanceChord :: Fractional a => [Waveform t a] -> Waveform t a
balanceChord notes = sampleFrom $ \t -> foldr (\x s -> s + factor * sample x t) 0 $ notes
  where
    factor = realToFrac . recip @Double . fromIntegral . length $ notes

-- | Play several waves on top of each other, without worrying about the volume. See 'balanceChord' for
-- a normalized version.
mergeWaves :: Num a => [Waveform t a] -> Waveform t a
mergeWaves notes = sampleFrom $ \t -> foldr (\x s -> s + sampleAt t x) 0 notes
  -- Average Frequency
  --,frequency = fmap (/(fromIntegral $ length notes)) . foldl (liftA2 (+)) (Just 0) . map frequency $ notes

-- | @'waveformToWAVE' outputLength@ gives a @'WAVE'@ file object by sampling the given @'DWave'@ at @44100Hz@.
-- May disbehave or clip based on behavior of @'doubleToSample'@ if the DWave takes values outside of @[-1,1]@.
waveformToWAVE :: Tick -> Int -> Wavetable -> WAVE.WAVE
waveformToWAVE outTicks sampleRate w = WAVE.WAVE
  {WAVE.waveHeader = WAVE.WAVEHeader
    {WAVE.waveNumChannels = 1
    ,WAVE.waveFrameRate = sampleRate
    ,WAVE.waveBitsPerSample = 32
    ,WAVE.waveFrames = Just $ fromIntegral outTicks
    }
  ,WAVE.waveSamples = [map (unDiscrete . sample w) [0 .. outTicks - 1]]
  }

-- | Take a specified number of samples at the given sample rate from the stream and put them in a @'WAVE.WAVE'@.
wavestreamToWAVE :: Tick -> Int -> [Discrete] -> WAVE.WAVE
wavestreamToWAVE outTicks sampleRate ws = WAVE.WAVE
  {WAVE.waveHeader = WAVE.WAVEHeader
    {WAVE.waveNumChannels = 1
    ,WAVE.waveFrameRate = sampleRate
    ,WAVE.waveBitsPerSample = 32
    ,WAVE.waveFrames = Just $ fromIntegral outTicks
    }
  ,WAVE.waveSamples = [take outTicks $ map unDiscrete ws]
  }

-- | Triangle wave of the given frequency
triWave :: Double -> Waveform Double Double
triWave f = sampleFrom $ \t -> let r = t * f - fromIntegral (floor (t * f) :: Int) in if r < 0.25
  then 4 * r
  else if r < 0.75
    then 2 - (4 * r)
    else -4 + (4 * r)

-- | Digital saw wave with given frequency
sawWave :: Double -> Waveform Double Double
sawWave f = sampleFrom $ \t -> 2 * (t * f - fromIntegral (floor (t * f) :: Int)) - 1

-- | Ramp from 0 to 1 over a time period
ramp :: Double -> Double -> Double
ramp = rampFrom 0

-- | Ramp from @x0@ to 1 from time zero to the max time
rampFrom :: Double -> Double -> (Double -> Double)
rampFrom x0 time t = if t < time
  then x0 + (1 - x0) * (max t 0)/time
  else 1

-- | Arbitrarily chosen standard tick rate, used in @'testWave'@
stdtr :: Num a => a
stdtr = 48000 --32000

-- | Output the first @len@ seconds of the given @'Wavetable'@ to a @.wav@ file at the given path for testing.
-- The volume is also attenuated by 50% to not blow out your eardrums.
-- Also pretty prints the wave.
testWave :: Double -> String -> Wavetable -> IO ()
testWave len fp w = print w >> pure w >>= WAVE.putWAVEFile (fp ++ ".wav") . waveformToWAVE (floor $ len*stdtr) stdtr . amplitudeModulate (sampleFrom $ const 0.5)

-- | Output the first @len@ seconds of the given wave stream to a @.wav@ file at the given path for testing.
testWaveStream :: Double -> String -> [Discrete] -> IO ()
testWaveStream len fp = WAVE.putWAVEFile (fp ++ ".wav") . wavestreamToWAVE (floor $ len*stdtr) stdtr

-- | Concert A4 frequency is 440Hz
concertA :: Num a => a
concertA = 440

-- | An envelope that suspends for @noteDuration@ before releasing
susEnvelope :: Envelope Tick Discrete -> Tick -> Wavetable
susEnvelope (Envelope del att hol dec sus rel) noteDuration = sampleFrom $ \t -> if t < del
  then 0
  else if t - del < att
    then doubleToDiscrete $ (t - del) `divTD` att
    else if t - del - att < hol
      then 1
      else if t - del - att - hol < dec
        then doubleToDiscrete $ 1 + ((t - del - att - hol) `divTD` dec) * (discreteToDouble sus - 1)
        else if t < noteDuration
          then sus
          else if t - noteDuration < rel
            then sus * doubleToDiscrete (1 - (t - noteDuration) `divTD` rel)
            else 0
  where
    divTD a b = fromIntegral a / fromIntegral b

-- | An envelope that suspends forever instead of releasing
suspendVelope :: Double -> Double -> Double -> Double -> Double -> DWave
suspendVelope del att hol dec sus = sampleFrom $ \t -> if t < del
  then 0
  else if t - del < att
    then (t - del) / att
    else if t - del - att < hol
      then 1
      else if t - del - att - hol < dec
        then 1 + (t - del - att - hol)/dec * (sus - 1)
        else sus

-- | Synthesize an envelope to a 'DWave'
envelope :: Envelope Double Double -> DWave
envelope (Envelope del att hol dec sus rel) = sampleFrom $ \t -> if t < del
  then 0
  else if t - del < att
    then (t - del) / att
    else if t - del - att < hol
      then 1
      else if t - del - att - hol < dec
        then 1 + (t - del - att - hol)/dec * (sus - 1)
        else if t - del - att - hol - dec < rel
          then sus * (1 - (t - del - att - hol - dec)/rel)
          else 0

-- | DAHDSR envelope with times @t@ and sustain level @b@
data Envelope t b = Envelope t t t t b t

-- | Discretize an envelope by converting times into ticks
discretizeEnvelope :: Double -> Envelope Double Double -> Envelope Tick Discrete
discretizeEnvelope tickRate (Envelope del att hol dec sus rel) = Envelope (dd del) (dd att) (dd hol) (dd dec) (doubleToDiscrete sus) (dd rel)
  where
    dd = floor . (*tickRate)

-- | Shift a wave in time to start at the specified time after its old start time
timeShift :: Num t => t -> Waveform t a -> Waveform t a
timeShift dt = sampleFrom . (. subtract dt) . sample

-- | Shift a wave in time such that the new zero is at the specified position
seekTo :: Num t => t -> Waveform t a -> Waveform t a
seekTo dt = sampleFrom . (. (+dt)) . sample

-- | Modify the amplitude of a wave by a constant multiple
setVolume :: Num a => a -> Waveform t a -> Waveform t a
setVolume = amplitudeModulate . sampleFrom . const

-- | The empty wave that is always zero when sampled
emptyWave :: Num a => Waveform t a
emptyWave = sampleFrom $ const 0

-- | Convolve with explicit discrete filter kernel weights.
discreteConvolve :: (Num a, Num t) => Waveform t [(t,a)] -> Waveform t a -> Waveform t a
discreteConvolve profile w = sampleFrom $ \t -> sum . map (\(dt,amp) -> amp * sample w (t + dt)) $ sample profile t

-- | This operation is not convolution, but something kind of like it. Use for creative purposes? Should be fast!
-- @wackyNotConvolution modf profile w = sampleFrom $ \t -> sample (modulate modf (sample profile t) w) t@
wackyNotConvolution :: (a -> b -> c) -> Waveform t (Waveform t a) -> Waveform t b -> Waveform t c
wackyNotConvolution modf profile w = sampleFrom $ \t -> sample (modulate modf (sample profile t) w) t

-- | Perform a discrete convolution. The output waveform is @f(t) = \int_{t-tickRadius}^{t+tickRadius} (kernel(t))(x) * w(t+x) dx@
-- but is discretized such that @x@ is always a multiple of @skipRate@.
tickConvolution :: (Show a,Fractional a)
                => Tick -- ^ @tickRadius@
                -> Tick -- ^ @skipRate@
                -> Waveform Tick (Waveform Tick a) -- ^ The kernel of the convolution at each @'Tick'@
                -> Waveform Tick a -- ^ w(t)
                -> Waveform Tick a
tickConvolution tickRadius skipRate profile w = sampleFrom $ \t -> let !kern = sample profile t in sum . map (\dt -> (*0.01) . (*sample w (t + dt)) . sample kern $ dt) $ sampleDeltas
  where
    sampleDeltas = map (*skipRate) [-stepsPerSide.. stepsPerSide]
    stepsPerSide = tickRadius `div` skipRate
    -- !stepModifier = realToFrac . recip . fromIntegral $ stepsPerSide :: Double

--accuCount :: IORef Integer
--accuCount = unsafePerformIO $ newIORef 0

-- | This computes a solid slice of a convolution much faster than sampling each point individually by using a convolution cache window.
fastTickConvolutionFixedKern :: Double -> Tick -> Tick -> Wavetable -> Wavetable -> [Discrete]
fastTickConvolutionFixedKern tickRate tickStart tickRadius kern w = kern `seq` go (vi,0) tickStart
  where
    vi = runST $ do
      v <- BA.newByteArray (windowSize * 4)
      forM_ [-tickRadius .. (tickRadius - 1)] $ \dt -> BA.writeByteArray v (1 + tickRadius + dt) (unDiscrete $ sample w (tickStart + dt))
      BA.unsafeFreezeByteArray v
    -- Table, breakpoint, offset
    stepSt :: (BA.ByteArray, Int) -> Tick -> (Discrete, BA.ByteArray)
    stepSt (!t,!b) o = runST $ do
      v <- BA.unsafeThawByteArray t
      BA.writeByteArray v (b `mod` windowSize) (unDiscrete $ sample w (o + tickRadius))
      let
        accu !tot (dt:dts) = do
          --let !() = unsafePerformIO $ readIORef accuCount >>= \x -> (if x `mod` 10000 == 0 then putStrLn ("accu: " ++ show x) else pure ()) >> writeIORef accuCount (x + 1)
          s <- BA.readByteArray v ((b + tickRadius + dt + 1) `mod` windowSize)
          let !toAdd = (0.01 * Discrete (BA.indexByteArray kernel dt)) * Discrete s
          accu (tot + toAdd) dts
        accu tot [] = pure tot
      k <- accu 0 sampleDeltas
      !v' <- BA.unsafeFreezeByteArray v
      pure (k,v')
    go (!t,!b) o = let (s,t') = stepSt (t,b) o in s `seq` (s : go (t',(b+1) `mod` windowSize) (o + 1))
    --accu :: BA.MutableByteArray s -> Int -> Tick -> Discrete -> Tick -> ST s Discrete
    windowSize = tickRadius * 2 + 1
    sampleDeltas = [-tickRadius .. tickRadius]
    !kernel = runST $ do
      v <- BA.newByteArray (windowSize * 4)
      forM_ [-tickRadius .. (tickRadius - 1)] $ \dt -> BA.writeByteArray v (1 + tickRadius + dt) (unDiscrete . (* doubleToDiscrete (recip tickRate)) $ sample kern dt)
      BA.unsafeFreezeByteArray v

-- | Same as @'tickConvolution'@ but for arbitarily valued waveforms. Works on @'DWave'@ for example.
sampledConvolution :: (RealFrac t, Fractional a)
                   => t -- ^ @convolutionSampleRate@, controls sampling for @x@
                   -> t -- ^ @convolutionRadius@, continuous analogue of @tickRadius@
                   -> Waveform t (Waveform t a) -- ^ Kernel of convolution for each time
                   -> Waveform t a -> Waveform t a
sampledConvolution convolutionSampleRate convolutionRadius profile w = sampleFrom $ \t -> sum . map (\dt -> (*(realToFrac . recip $ convolutionSampleRate * convolutionRadius)) . (* sample w (t + dt)) . sample (sample profile t) $ dt) $ sampleDeltas
  where
    sampleDeltas = map ((/convolutionSampleRate) . realToFrac) [-samplesPerSide .. samplesPerSide]
    samplesPerSide = floor (convolutionRadius * convolutionSampleRate) :: Int
    --sampleCount = 2 * samplesPerSide + 1

-- | Makes a filter which selects frequencies near @bandCenter@ with tuning parameter @bandSize@.
-- Try: @'optimizeFilter' 200 . 'tickTable' 'stdtr' $ 'bandpassFilter' 'concertA' 100@
bandpassFilter :: Fractional a
               => Double -- ^ @bandCenter@
               -> Double -- ^ @bandSize@
               -> Waveform Double a
bandpassFilter bandCenter bandSize = sampleFrom $ \t -> if t == 0 then 1 else realToFrac $ (sin (bandFreq * t)) / (bandFreq * t) * (cos (centerFreq * t))
  where
    !bandFreq = 2 * pi * bandSize
    !centerFreq = 2 * pi * bandCenter

-- | Ladder filter, adapted from VCV Rack
ladderFilter :: Double -> [Double] -> [Double] -> [Double] -> [(Double,Double,Double,Double)]
ladderFilter sampleRate reso cutoffs inputs = fmap (FV.convert :: FV.VecList 4 Double -> (Double,Double,Double,Double)) $ rkSolveStepSize (1/sampleRate) rkf (FV.mk4 1e-6 1e-6 1e-6 1e-6) (zip3 reso cutoffs inputs)
  where
    rkf (r,c,i) x = let cut = 2 * pi * c in FV.mk4 (cut * (i - r * (x FV.! 2) - (x FV.! 0))) (cut * (x FV.! 0 - x FV.! 1)) (cut * (x FV.! 1 - x FV.! 2)) (cut * (x FV.! 2 - x FV.! 3))

-- | Use Runge-Kutta method to solve the diffeq
rkSolve
  :: Applicative f => (i -> f Double -> f Double) -- ^ User suplied function mapping inputs to a coefficient matrix
  -> f Double -- ^ Initial value
  -> [i] -- ^ Stream of inputs
  -> [f Double] -- ^ Stream of outputs
rkSolve rkf = scanl stepSt
  where
    stepSt x i = let
      k1 = rkf i x
      k2 = rkf i (liftA2 (\k xi -> xi + k/2) k1 x)
      k3 = rkf i (liftA2 (\k xi -> xi + k/2) k2 x)
      k4 = rkf i (liftA2 (+) k3 x)
      in liftA2 (\xi ki -> xi + ki/6) x $ liftA2 (+) (liftA2 (+) (fmap (*2) k2) (fmap (*2) k3)) (liftA2 (+) k1 k4)

-- | Like @'rkSolve'@ but using non-unit step size by multiplying the coefficient matrix by a factor
rkSolveStepSize :: Applicative f => Double -> (i -> f Double -> f Double) -> f Double -> [i] -> [f Double]
rkSolveStepSize h rkf = rkSolve ((fmap (* h) .) . rkf)

-- | Drive scaling function adapted from VCV Rack
scaledDrive :: Double -> Double
scaledDrive drv = (1 + drv) ** 5

-- | Uses @'ladderFilter'@ to low pass filter the input with the given resonance and cutoff
lowPassFilter :: Double -> [Double] -> [Double] -> [Double] -> [Double]
lowPassFilter rate res cut = fmap (\(_,_,_,d) -> d) . ladderFilter rate res cut

-- | Synthesize from a continuous profile by sampling frequencies.
-- Doesn't work very well; prefer 'synthFromDisreteProfile'.
synthFromFreqProfile :: (Double,Double) -> Double -> DWave -> DWave
synthFromFreqProfile (f0,f1) fSampRate prof = sampleFrom $ \t -> max (-1) $ min 1 $ sum . fmap ((/fromIntegral (length sampPoints)) . sampleAt t . getComponent) $ sampPoints
  where
    getComponent f = fmap (*sample prof f) $ sinWaveWithPhase f (pi/2)
    sampPoints = [f0, f0 + 1/fSampRate .. f1]
    --jitterFactor k = (fromIntegral (floor (((k - f0) / (f1 - f0)) * 81083) `mod` 115)) / 115
    --jitter k = k + (2/fSampRate) * jitterFactor k

-- | Continuous saxopohone profile adapted from https://asa.scitation.org/doi/10.1121/1.396474
saxProfile :: DWave
saxProfile = sampleFrom $ \f -> let x = f/fb in sqrt (n * x/(1 + x**7))
  where
    fb = 618
    n = 1.088910

-- | Synthesize a wavetable from the specified relative strengths of various frequencies.
-- Essentially an inverse fourier transform
synthFromDiscreteProfile :: [(Double,Double)] -> Wavetable
synthFromDiscreteProfile fs = {-solidSlice 0 (truncate $ stdtr/(1 :: Double)) .-} tickTable stdtr . discretize . mergeWaves . fmap (\(f,amp) -> fmap (*(normFactor * amp)) $ sinWaveWithPhase f f) $ fs
  where
    normFactor = 0.9/(sum $ fmap snd fs)

-- | Produces a discrete profile of relative harmonic strengths from the given functions.
-- Try @'harmonicEquationToDiscreteProfile' (\x -> 0.1894 / (x ** 1.02)) (\x -> 0.0321 / (x ** 0.5669))@
harmonicEquationToDiscreteProfile :: (Double -> Double) -> (Double -> Double) -> Double -> [(Double,Double)]
harmonicEquationToDiscreteProfile oddF evenF f0 = go (15 :: Int)
  where
    go 0 = []
    go k = (fromIntegral k * f0,(if even k then evenF else oddF) $ fromIntegral k) : go (k-1)

-- | Roughly resemebles the spectrum of a saxophone
discSaxProfile :: [(Double,Double)]
discSaxProfile = f1 ++ f2
  where
  f1 =
    [(588.6,0.0911)
    ,(1179,0.05463)
    ,(1767,0.03471)
    ,(2357,0.01109)
    ,(2946,0.01625)
    ] -- y=63/x ish
  f2 =
    [(294.9,0.03032)
    ,(884,0.02445)
    ,(1474,0.01172)
    ,(2062,0.008224)
    ,(2651,0.00881)
    ,(3241,0.005567)
    ]

-- | Roughly resemebles the spectrum of a saxophone
genSaxProfile' :: Double -> [(Double,Double)]
genSaxProfile' f0 =
  [(1 * f0,0.0474475510012406)
  ,(2 * f0,0.0141427241744339)
  ,(3 * f0,0.00838229600425522)
  ,(4 * f0,0.00452452580892121)
  ,(5 * f0,0.0129662632033567)
  ,(6 * f0,0.0102954061023286)
  ,(7 * f0,0.00246665413386754)
  ,(8 * f0,0.00309631105551752)
  ,(9 * f0,0.00313836031805197)
  ,(10 * f0,0.00549688733850934)
  ,(11 * f0,0.00826471159302181)
  ,(12 * f0,0.00280485189554740)
  ]

-- | Roughly resemebles the spectrum of a saxophone
genSaxProfile :: Double -> [(Double,Double)]
genSaxProfile f0 = f1 ++ f2
  where
  f1 =
    [(2 * f0,0.0911)
    ,(4 * f0,0.05463)
    ,(4 * f0,0.03471)
    ,(6 * f0,0.01109)
    ,(8 * f0,0.01625)
    ] -- y=63/x ish
  f2 =
    [(f0,0.03032)
    ,(3 * f0,0.02445)
    ,(5 * f0,0.01172)
    ,(7 * f0,0.008224)
    ,(9 * f0,0.00881)
    ,(11 * f0,0.005567)
    ]

-- | Rough saxophone timbre from @'gen
fakedSaxTimbre :: Double -> Wavetable
fakedSaxTimbre = synthFromDiscreteProfile . genSaxProfile'

-- | A timbre that consists of several sin wave voices at fixed chord intervals from each other
chordSinTimbre :: ChordVoicing PitchFactorDiagram -> Double -> Wavetable
chordSinTimbre c r = discretize . tickTable stdtr . balanceChord . fmap (sinWave . flip intervalOf r) $ getVoiceList c

-- | Discretize the output of a @'Double'@ producing waveform
discretize :: Functor f => f Double -> f Discrete
discretize = fmap doubleToDiscrete

-- | Discretize the input to a @'Double'@ consuming waveform
tickTable :: Double -- ^ Sample rate. Each tick is @1/sampleRate@ seconds
          -> Waveform Double a -> Waveform Tick a
tickTable tickrate w = sampleFrom $ \t -> sample w (fromIntegral t/tickrate)

-- | Discretize the input to a @'Double'@ consuming waveform
tickTablePer :: Double -- ^ Sample period. Sampling rate is @1/period@
          -> Waveform Double a -> Waveform Tick a
tickTablePer per w = sampleFrom $ \t -> sample w (fromIntegral t * per)

-- | A domain- and codomain-discretized @'Waveform'@ suitable for writing to a WAVE file.
-- See @'waveformToWAVE'@.
type Wavetable = Waveform Tick Discrete

-- | A data structure for storing the results of a @'Wavetable'@ on some subset of its domain.
-- Used internally.
data CompactWavetable = CompactWavetable {getWavetable :: Vector.Vector Int32}

fromCompact :: CompactWavetable -> Wavetable
fromCompact cwt = sampleFrom $ \t -> case getWavetable cwt Vector.!? t of
  Just d -> Discrete d
  Nothing -> error "fromCompact Wavetable sampled outside size of compact table"

-- | Optimize a @'Wavetable'@ by storing its values in a particular range.
-- Uses @(tickEnd - tickStart + 1) * sizeOf (_ :: 'Discrete')@ bytes of memory to do this.
solidSlice :: Tick -> Tick -> Wavetable -> Wavetable
solidSlice tickStart tickEnd w = sampleFrom $ \t -> case getWavetable cwt Vector.!? (t-tickStart) of
  Just d -> Discrete d
  Nothing -> sample w t
  where
    cwt = CompactWavetable {getWavetable = Vector.generate (tickEnd - tickStart + 1) (unDiscrete . sample w . (+tickStart))}

-- | Optimize a filter by doing @'solidSlice'@ around @t=0@ since those values are sampled repeatedly in a filter
optimizeFilter :: Tick -> Wavetable -> Wavetable
optimizeFilter tickRadius = solidSlice (-tickRadius) tickRadius

-- | Take the Fourier Transform of a complex valued @'Tick'@ sampled waveform
fourierTransform :: Tick -> Double -> Waveform Tick (Complex Double) -> Waveform Double (Complex Double)
fourierTransform tickRadius fTickRate x = sampleFrom $ \f -> sum . map (\n -> sample x n / (fromIntegral tickRadius) * cis (2 * pi * f * (fromIntegral n / fTickRate))) $ [-tickRadius .. tickRadius]

-- | Take the Fourier Transform of a @'Wavetable'@
realDFT :: Tick -- ^ Radius of Fourier Transform window in @'Tick'@s. Try 200
        -> Double -- ^ Sampling rate to use for the Fourier transform. Try the sample sample rate as the @'Wavetable'@
        -> Wavetable -> Wavetable
realDFT tickRadius fTickRate w = discretize $ tickTable 1 $ fmap ((min 1) . magnitude) $ fourierTransform tickRadius fTickRate ((\x -> discreteToDouble x :+ 0) <$> solidSlice (-tickRadius) tickRadius w)

-- | Skip every @n@ ticks in the in the given @'Waveform'@.
-- @'sample' ('skipTicks' n w) k = 'sample' w (n*k)@
skipTicks :: Tick -- ^ @n@
          -> Waveform Tick a -> Waveform Tick a
skipTicks skipRate w = sampleFrom $ \t -> sample w (skipRate * t)

skipStream :: Tick -> [Discrete] -> [Discrete]
skipStream skipRate (x:xs) = x : skipStream skipRate (drop skipRate xs)
skipStream _ [] = []

-- | Optimize a @'Wavetable'@ that we know to be periodic by storing it's values on one period.
-- Takes @period * sizeOf (_ :: 'Discrete')@ bytes of memory to do this.
--
-- Warning: Causes clicking on period boundaries if the period isn't exactly the given value in ticks.
exploitPeriodicity :: Tick -- ^ Period in @'Tick'@s of the @'Wavetable'@.
                   -> Wavetable -> Wavetable
exploitPeriodicity period x = sampleFrom $ \t -> case getWavetable cwt Vector.!? (t `mod` period) of
  Just d -> Discrete d
  Nothing -> sample x t
  where
    cwt = CompactWavetable {getWavetable = Vector.generate period (unDiscrete . sample x)}

-- | Attempts to do a fast fourier transform, but the units of the domain of the output are highly suspect.
-- May be unreliable, use with caution.
usingFFT :: Tick -> Wavetable -> Wavetable
usingFFT tickRadius w = sampleFrom $ \t -> if t < (fromIntegral $ length l)
  then (!! t) . fmap (doubleToDiscrete . magnitude) $ l
  else 0
  where
    l = fft (map ((\x -> discreteToDouble x :+ 0) . sample w) [-tickRadius .. tickRadius])

-- | Cooley-Tukey fft
fft :: [Complex Double] -> [Complex Double]
fft [] = []
fft [x] = [x]
fft xs = zipWith (+) ys ts ++ zipWith (-) ys ts
    where n = length xs
          ys = fft evens
          zs = fft odds
          (evens, odds) = split xs
          split [] = ([], [])
          split [x] = ([x], [])
          split (x:y:xss) = (x:xt, y:yt) where (xt, yt) = split xss
          ts = zipWith (\z k -> exp' k n * z) zs [0..]
          exp' :: Int -> Int -> Complex Double
          exp' k a = cis $ -2 * pi * (fromIntegral k) / (fromIntegral a)

max_stars_repo_name: Lazersmoke/boopadoop

max_issues_repo_name: Lazersmoke/boopadoop

max_forks_repo_name: Lazersmoke/boopadoop


text:
{-# LANGUAGE BangPatterns, DeriveFunctor, DeriveFoldable, DeriveTraversable, DeriveDataTypeable #-}
-- | A scene graph is used to represent any N-dimensional scene. Since it is
--   polymorphic over its contents, we can use it to represent a variety of
--   things, from a physics world to a render tree.
module Game.SceneGraph ( -- * Basic Data Types
                         SceneGraph(..)
                       , TransformView(..)
                       -- * Views
                       , viewTrans
                       , matRep
                       , Quaternion
                       -- * Transformations
                       , translate
                       , scaleGraph
                       , rotate2D
                       , rotate3D
                       , rotateQuat
                       , rotate3DX
                       , rotate3DY
                       , rotate3DZ
                       , rotateEuler
                       , applyPoint
                       , applyVector
                       -- * Graph Visitors
                       , optimize
                       , sceneMap
                       ) where

import Prelewd hiding (toList, join, (<>))

import Impure

import Control.DeepSeq
import Control.Monad.Par
import Data.Function.Pointless
import Data.Typeable
import Storage.List (init)
import Numeric.LinearAlgebra
import Text.Show

-- | Represents an N-dimensional scene graph, with objects at the leaves, and
--   transformations applied in a hierarchical fashion.
--
--   This is an instance of Functor, Foldable, and Traversable, so feel free to
--   walk it like you would any recursive data structure.
data SceneGraph a = Object a
                  | Transform !Transformation !(SceneGraph a)
                  | Branch [SceneGraph a]
    deriving (Show, Functor, Foldable, Typeable)

instance NFData a => NFData (SceneGraph a) where
    rnf (Transform t g) = rnf t `seq` rnf g
    rnf (Object x)      = rnf x
    rnf (Branch xs)     = rnf xs

-- | We represent the basic transformations (translate, scale, rotate) as their
--   "easy" representation, as well as the more complicated matrix. This allows
--   us to optimize the scene graph with the semantic information of chained
--   transformations instead of just blindly multiplying matricies.
--
--   For example, instead of multiplying two (rather complicated) quaternion
--   matricies to chain rotation, we just perform quaternion multiplication.
--
--   By keeping the matricies around too, we can always fall back on matrix
--   multiplication, and we can also avoid recomputing them.
--
--   Oh yeah, and all the matricies use homogenous coordinates, so they're
--   (N+1)x(N+1) in an n-dimensional scene. The more you know.
data Transformation = Translate !(Vector Double)
                                 (Matrix Double)
                    | Scale     !(Vector Double)
                                 (Matrix Double)
                    | Rotate2D  !Double -- θ - radians, clockwise.
                                 (Matrix Double) -- affine representation
                    | AffineTransform !(Matrix Double) -- fall back generic transformations
                                                       -- also used for rotations.
    deriving (Show, Typeable)

-- | A data type used for deconstructing transformations and seeing what they
--   "really are".
data TransformView = TranslateView !(Vector Double)
                   | ScaleView !(Vector Double)
                   | Rotate2DView !Double -- ^ θ - radians, clockwise.
                   | AffineTransformView !(Matrix Double) --  ^ Uses homogenous coordinates.

-- | Use this function to deconstruct a transformation into what it represents.
--
--   We use a more complicated internal data type for representing
--   transformations, and this view can be used to simplify it. The GHC extension
--   "ViewPatterns" will probably be of help here.
--
--   > -- f returns a scaled (by a factor of 2) translation vector.
--   > f :: Transformation -> Vector Double
--   > f (viewTrans -> Translate v) = v*2
viewTrans :: Transformation -> TransformView
viewTrans (Translate v _)     = TranslateView v
viewTrans (Scale v _)         = ScaleView v
viewTrans (Rotate2D θ _)      = Rotate2DView θ
viewTrans (AffineTransform m) = AffineTransformView m

-- | Retrieves the matrix view of a transformation. This will be an (N+1)x(N+1)
--   matrix, where N is the dimensions in the scene.
matRep :: Transformation -> Matrix Double
matRep (Translate _ m)     = m
matRep (Scale _ m)         = m
matRep (Rotate2D _ m)      = m
matRep (AffineTransform m) = m

seq_ :: a -> ()
seq_ = (`seq` ())

instance NFData Transformation where
    rnf (Translate _ m)     = seq_ m
    rnf (Scale     _ m)     = seq_ m
    rnf (Rotate2D  _ m)     = seq_ m
    rnf (AffineTransform _) = ()

-- | Just a 4-element vector, specifying a 3-dimensional rotation. See the
--   wikipedia page for more information on Quaternions.
type Quaternion = Vector Double

-- | Represents a 2-dimensional, clockwise rotation by θ radians. This is
--   probably easier to use than building a rotation matrix manually.
rotate2D :: Double -- ^ θ.
         -> Transformation
rotate2D !θ = Rotate2D θ $ rotate2DMat θ

rotate2DMat :: Double
            -> Matrix Double
rotate2DMat !θ = (3><3) [  c, s, 0
                        , -s, c, 0
                        , 0,  0, 1 ]
    where
        c = cos θ
        s = sin θ

-- | Rotates by a given quaternion.
--
--   Note: This algorithm is untested, and disagrees with wikipedia right now.
--         I'll look into it one day.
rotateQuat :: Quaternion -> Transformation
rotateQuat q = AffineTransform $ (4><4) [ 1-2.0*(yy-zz), 2.0*(xy-zw)    , 2.0*(xz+yw), 0
                                        , 2.0*(xy+zw)  , 1.0-2.0*(xx-zz), 2.0*(yz-xw), 0
                                        , 2.0*(xz-yw)  , 2*(xw+yz)      , 1-2*(xx-yy), 0
                                        , 0            , 0              , 0          , 1 ]
    where
        [ w, x, y, z ] = toList q
        xx = x*x
        yy = y*y
        zz = z*z
        zw = z*w
        xy = x*y
        xz = x*z
        xw = x*w
        yw = y*w
        yz = y*z

-- | Rotates by a given angle around an arbitrary 3-dimensional axis.
rotate3D :: Double -- ^ The amount (in radians) to rotate by.
         -> Vector Double -- ^ The 3-dimensional axis around which we are rotating.
         -> Transformation
-- https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
rotate3D !θ !axis = AffineTransform $ (4><4) [ c+x*xoc , x*yoc-zs, x*zoc+ys, 0
                                             , y*xoc+zs, c+y*yoc , y*zoc-xs, 0
                                             , z*xoc-ys, z*yoc+xs, c+z*zoc , 0
                                             ,     0   ,    0    ,    0    , 1 ]
    where
        c = cos θ
        oc = 1-c
        s = sin θ
        [x, y, z] = toList axis
        xs = x*s
        ys = y*s
        zs = z*s
        xoc = x*oc
        yoc = y*oc
        zoc = z*oc

-- | Rotates θ radians clockwise around the X-axis.
rotate3DX :: Double -- θ
          -> Transformation
rotate3DX !θ = AffineTransform $ (4><4) [ 1,  0, 0, 0
                                        , 0,  c, s, 0
                                        , 0, -s, c, 0
                                        , 0,  0, 0, 1 ]
    where
        c = cos θ
        s = sin θ

-- | Rotates θ radians clockwise around the Y axis.
rotate3DY :: Double -- θ
          -> Transformation
rotate3DY !θ = AffineTransform $ (4><4) [  c, 0, -s, 0
                                        ,  0, 1,  0, 0
                                        , -s, 0,  c, 0
                                        ,  0, 0,  0, 1 ]
    where
        c = cos θ
        s = sin θ

-- | Rotates θ radians clockwise around the Y axis.
rotate3DZ :: Double -- θ
          -> Transformation
rotate3DZ !θ = AffineTransform $ (4><4) [  c,  s, 0, 0
                                        , -s,  c, 0, 0
                                        ,  0,  0, 1, 0
                                        ,  0,  0, 0, 1 ]
    where
        c = cos θ
        s = sin θ

-- | Rotates around a yaw, pitch, and roll - all in one transformation.
rotateEuler :: Double -- X-axis rotation.
            -> Double -- Y-axis rotation.
            -> Double -- Z-axis rotation.
            -> Transformation
rotateEuler x y z = let (AffineTransform x') = rotate3DX x
                        (AffineTransform y') = rotate3DY y
                        (AffineTransform z') = rotate3DZ z
                     in AffineTransform $ z' * y' * x'

-- | Adds a number to the end of a vector. Tends to be useful when constructing
--   homogenous coordinates.
vappend :: Element a => a -> Vector a -> Vector a
vappend x = join . flip (:) [ constant x 1 ]

-- | Applies the given modelview matrix to the given vector.
applyVector :: Matrix Double -> Vector Double -> Vector Double
applyVector m v = subVector 0 (dim v) $ m <> vappend 1.0 v

-- | Applied the given modelview matrix to the given point.
applyPoint :: Matrix Double -> Vector Double -> Vector Double
applyPoint m p = subVector 0 (dim p) $ m <> vappend 0.0 p

scaleMat :: Vector Double -> Matrix Double
scaleMat = diag . vappend 1.0

-- | Performs an n-dimensional scaling of the scene graph. Bigger numbers mean
--   bigger objects. 1.0 means no scaling.
scaleGraph :: Vector Double -- ^ The N-dimensional vector of values to scale by.
           -> Transformation
scaleGraph v = Scale v $ scaleMat v

translateMat :: Vector Double -> Matrix Double
translateMat v = fromColumns $ ((<?> error "Empty vector") . init . toColumns . ident $ dim v + 1) `mappend` [vappend 1.0 v]

-- | Translates an n-dimensional scene graph by the given amount. Bigger values
--   mean bigger objects. 1.0 means no scaling.
translate :: Vector Double -- ^ The N-dimensional vector of values to translate by.
          -> Transformation
translate v = Translate v $ translateMat v

-- | Cleans up a scene graph, and does its best to avoid as many
--   multiplications and indirections as possible.
--
--   This will walk the whole graph, but if you're doing so anyhow, it'd be a
--   great way to avoid a lot of extra processing due to inefficient structure.
optimize :: SceneGraph a -> SceneGraph a
optimize x@(Object _) = x
optimize (Branch xs)  = optBranch xs
    where
        optBranch :: [SceneGraph a] -> SceneGraph a
        optBranch []  = Branch [] -- lolwut
        optBranch [x] = optimize x
        optBranch xs' = Branch $ foldr f [] xs'
            where
                -- merges sub-branches up into their parent.
                f y ys = case optimize y of
                            Branch zs -> zs `mappend` ys
                            g         -> g:ys
optimize (Transform t g) = go
    where
        go = let subGraph = optimize g
              in case subGraph of
                Transform t' subsubGraph -> case (t, t') of
                    (Translate v _, Translate v' _) -> Transform (combineTrans v v') subsubGraph
                    (Scale v _, Scale v' _)         -> Transform (combineScale v v') subsubGraph
                    (Rotate2D θ _, Rotate2D τ _)    -> Transform (combineRot θ τ) subsubGraph
                    (_, _)                          -> Transform t subGraph
                _ -> Transform t subGraph
        combineTrans = translate .: (+)
        combineScale = scaleGraph .: (+)
        combineRot   = rotate2D .: (+)

-- | Maps a function taking an object and its modelview matrix, returning a new
--   scene graph of a possibly different type.
--
--   Please optimize the scene graph before calling this function for optimal
--   efficiency.
--
--   To transform a point by the given matrix, use 'homPoint'. To
--   transform a vector, use 'homVector'.
--
--   Don't do the multiplication yourself. You'll probably get it wrong.
sceneMap :: NFData b
         => Int -- ^ The number of dimensions in the scene.
         -> (a -> Matrix Double -> b) -> SceneGraph a -> Par (SceneGraph b)
sceneMap n f = sMap (ident n)
    where
        -- Give a stack of matricies to multiply and a scene graph, returns the
        -- transformed (by f) scene graph.
        sMap m (Object x)  = return . Object $ f x m
        sMap m (Transform t child) = Transform t <$> sMap (matRep t * m) child
        sMap m (Branch xs) = Branch <$> parMapM (sMap m) xs
{-# INLINE sceneMap #-}

max_stars_repo_name: bfops/Chess

max_issues_repo_name: bfops/Chess

max_forks_repo_name: bfops/Chess


text:
{-# language TypeApplications #-}

module BaseSpec where

import Data.Complex
import Data.Fixed (mod')
import Data.Semigroup
import Geometry
import Geometry.Testing
import Test.Hspec
import Test.Invariant
import Test.QuickCheck hiding (scale)

spec :: Spec
spec = do
  describe "Direction" $ do
    describe "equality" $ do
      it "1" $ asDeg 0 == asDeg 360
      it "2" $ asDeg 10 == asDeg 370
      it "3" $ 45 == asCmp @Direction (1 :+ 1)
      it "4" $ 90 == asCmp @Direction (0 :+ 1)
      it "5" $ 90 == asDeg (90 + 1e-12)
      it "6" $ asCmp @Direction (2 :+ 3) == asCmp @Direction (4 :+ 6)

    describe "inequality" $ do
      it "1" $ asDeg 0 <= asDeg 360
      it "2" $ asDeg 0 < asDeg 360.0001
      it "3" $ property $ \a -> 0 <= turns a && turns a <= 1

    describe "isomorphism" $ do
      it "asRad" $ property $ \x -> (rad . asRad) x ~= x `mod'` (2*pi)
      it "asDeg" $ property $ \x -> (deg . asDeg) x ~= x `mod'` 360
      it "asTurns" $ property $ \x -> (turns . asTurns) x ~= x `mod'` 1

    describe "radians" $ do
      it "1" $ rad 10 ~= rad 370
      it "2" $ rad (-10) ~= rad 350

  describe "Affine" $ do
    describe "isomorphism" $ do
      it "1" $ property $ cmp `inverts` xy
      it "2" $ property $ xy `inverts` cmp

    describe "dot" $ do
      it "pair" $ dot @XY @XY (1, 2) (3, 6) == 15
      it "complex" $ dot @Cmp @Cmp (1 :+ 2) (3 :+ 6) == 15
      it "orthogonality" $ property $ \x ->
        not (isZero x) ==> rotate 90 x `isOrthogonal` (x :: XY)

    describe "cross" $ do
      it "pair" $ cross @XY @XY (1, 2) (3, 6) == 0
      it "complex" $ cross @Cmp @Cmp (1 :+ 2) (3 :+ 6) == 0

    describe "norm" $ do
      it "pair" $ norm @XY (3, 4) == 5
      it "complex" $ norm @Cmp (3 :+ 4) == 5

    describe "normalize" $ do
      it "pair" $ norm (normalize @XY (3, 4)) == 1
      it "complex" $ norm (normalize @Cmp (3 :+ 4)) == 1
      it "trivial" $ normalize @Cmp 0 == 0

    describe "opposite" $
      it "1" $ property $ \x -> x `isOpposite` opposite @XY x

  describe "Transformations" $ do
    describe "translate" $ do
      it "1" $ translate' @XY @XY (10, 20) (3, 4)  == (13, 24)
      it "2" $ translate' @XY @Cmp 10 (3, 4)  == (13, 4)
      it "3" $ translate' @Cmp @XY (3, 4) 10  == 13 :+ 4

    it "scale" $ scale @XY 10 (3, 4)  == (30, 40)
    it "scaleX" $ scaleX @XY 10 (3, 4) ~= (30, 4)
    it "scaleY" $ scaleY @XY 10 (3, 4) ~= (3, 40)

    describe "scaleAt" $ do
      it "1" $ scaleAt' @XY @XY (3, 4) 10 (3, 4)  == (3, 4)
      it "2" $ scaleAt' @XY @XY (0, 0) 10 (3, 4)  == (30, 40)
      it "3" $ scaleAt' @XY @XY (1, 1) 10 (3, 4)  == (21, 31)
      it "4" $ scaleAt' @XY @XY (1, 1) 10 (0, 0)  == (-9, -9)

    describe "rotate" $ do
      it "1" $ rotate @XY 0 (3, 4)  == (3, 4)
      it "2" $ rotate @XY 90 (3, 4)  ~= (-4, 3)
      it "3" $ rotate @XY 360 (3, 4)  == (3, 4)
      it "4" $ rotate @XY (asCmp (0 :+ 1)) (3, 4) ~= (-4, 3)

    describe "rotateAt" $ do
      it "1" $ rotateAt' @XY @XY (0, 0) 0 (3, 4) == (3, 4)
      it "2" $ rotateAt' @XY @XY (3, 4) 90 (3, 4) == (3, 4)
      it "3" $ rotateAt' @XY @XY (1, 0) 90 (0, 0) ~= (1, -1)
      it "4" $ rotateAt' @XY @XY (1, 1) 180 (0, 0) ~= (2, 2)

  describe "distanceTo" $ do
    let c1 = aCircle
        c2 = aCircle # scale 2
        c3 = aCircle # scale 3
    it "1" $ property $ \x -> distanceTo (c1 @-> x) c2 ~= 1
    it "2" $ property $ \x -> distanceTo (c3 @-> x) c2 ~= 1
    it "3" $ distanceTo aPoint c2 == 2
    it "4" $ distanceTo aPoint c3 == 3
    it "5" $ property $ \x a ->
      let l1 = aLine # rotate a
          l2 = l1 # translate' (normal l1 0)
      in distanceTo (l1 @-> x) l2 ~= 1
    it "6" $ property $ \x (Positive n) ->
      let t1 = regularPoly (n + 3)
          t2 = t1 # scale 2
      in distanceTo (t1 @-> x) t2 ~= cos (pi/ fromIntegral (n + 3))

  describe "Box" $ do
    let l = Segment (1:+2, 4:+6)
    it "0" $ box l == ((1,2),(4,6))
    it "1" $ (left.lower.corner $ l) == 1:+2
    it "2" $ (right.lower.corner $ l) == 4:+2
    it "3" $ (left.upper.corner $ l) == 1:+6
    it "4" $ (right.upper.corner $ l) == 4:+6
    it "5" $ figureHeight l == 4
    it "6" $ figureWidth l == 3

  describe "AlmostEqual" $ do
    it "1" $ property $ \i j -> let _ = (i :: Int, j :: Int)
                                in (i == j) == (i ~= j)
    it "2" $ property $ \i j -> let _ = (i :: Maybe Int, j :: Maybe Int)
                                in (i == j) == (i ~= j)

max_stars_repo_name: MethaHardworker/geometry

max_issues_repo_name: MethaHardworker/geometry

max_forks_repo_name: MethaHardworker/geometry