text:
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE DataKinds         #-}
{-# LANGUAGE TypeFamilies      #-}
{-# LANGUAGE RecordWildCards   #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GADTs             #-}
{-# LANGUAGE FlexibleContexts  #-}

module Main where

import           Lens.Micro.Extras (view)
import qualified Data.Text.IO as Text (readFile)
import           Data.Text (lines)
import           Data.Maybe (mapMaybe, fromMaybe)
import           Data.List (zip4) --, foldl')

import           Control.Monad (forM_)
import           Control.Monad.Random (evalRand)
import           System.Random (mkStdGen, RandomGen)
import qualified System.Random as SR (split)

import           Data.Semigroup ((<>))
import           Options.Applicative (Parser, option, auto, optional, long, short, value,
                                      strOption, helper, idm, execParser, info, (<**>), (<|>),
                                      flag')
import qualified Data.ByteString as B
import           Data.Serialize (Serialize, runPut, put, runGet, get) --, Get
import           System.Directory (createDirectoryIfMissing)

import           Numeric.LinearAlgebra.Static ((#), (&), R, unwrap, konst) --, zipWithVector)
import           Numeric.LinearAlgebra (toList)

--import           GHC.TypeLits
import           Data.Singletons (Sing, SomeSing(SomeSing), toSing)
import           Data.Singletons.TypeLits (KnownNat, Nat, Sing(SNat))
import           Data.Singletons.Prelude.List (Head, Last)

import           Grenade (Network,  FullyConnected, Tanh, Logit, Shape(D1), --Relu,
                         LearningParameters(LearningParameters), randomNetwork, S(S1D),
                         Gradients)
--import           Grenade (Network((:~>), NNil),  FullyConnected(FullyConnected), FullyConnected'(FullyConnected'), Tanh, Logit, Shape(D1), runNet, Relu,
--                         LearningParameters(LearningParameters), randomNetwork, S(S1D))

import           GrenadeExtras (binaryNetError, trainOnBatchesEpochs, normalize, hotvector)
import           GrenadeExtras.Zip (Zip)
--import           GrenadeExtras.Orphan()
import           GrenadeExtras.GradNorm (GradNorm)

import           Song (Song, line2song, genre)
import           Song.Grenade (song2TupleRn) --, SongSD)
import           Shuffle (shuffle) -- TODO: use the "standard" shuffleM http://hackage.haskell.org/package/random-shuffle

--import           Debug.Trace (trace)

type LogisticRegression = Network '[ FullyConnected 75 1, Logit ]
                                  '[ 'D1 75, 'D1 1, 'D1 1 ]

type OneHiddenLayer n = Network '[ FullyConnected 75 n, Tanh, FullyConnected n 1, Logit ]
                                '[ 'D1 75, 'D1 n, 'D1 n, 'D1 1, 'D1 1 ]

type DiscreteSepContinuousNet = Network '[ Zip ('D1 21) ('D1 2) (FullyConnected 21 2) ('D1 54) ('D1 2) (FullyConnected 54 2),
                                           Tanh, FullyConnected 4 1,
                                           Logit ]
                                        '[ 'D1 75, 'D1 4, 'D1 4, 'D1 1, 'D1 1 ]

type HugeNetwork = Network '[ FullyConnected 75 40, Tanh, FullyConnected 40 10, Tanh, FullyConnected 10 3, Tanh, FullyConnected 3 1, Logit ]
                           '[ 'D1 75, 'D1 40, 'D1 40, 'D1 10, 'D1 10, 'D1 3, 'D1 3, 'D1 1, 'D1 1 ]

--type FFNet = Network '[ FullyConnected 75 100, Tanh, FullyConnected 100 40, Tanh, FullyConnected 40 20, Relu, FullyConnected 20 1, Logit ]
--                     '[ 'D1 75, 'D1 100, 'D1 100, 'D1 40, 'D1 40, 'D1 20, 'D1 20, 'D1 1, 'D1 1 ]

--netLoad :: KnownNat n => FilePath -> IO (OneHiddenLayer n)
netLoad :: Serialize b => FilePath -> IO b
netLoad modelPath = do
  modelData <- B.readFile modelPath
  either fail return $ runGet get modelData

createKFoldPartitions :: Int -> [a] -> (Int, [([a],[a])])
createKFoldPartitions size ks = (parts, take parts . separateTrainAndValidLists . splitInSizes $ ks)
  where
    len = length ks
    parts = floor $ fromIntegral len / (fromIntegral size :: Double)

    splitInSizes :: [a] -> [[a]]
    splitInSizes [] = []
    splitInSizes xs =
      let (start, finish) = splitAt size xs
       in start : splitInSizes finish

    separateTrainAndValidLists :: [[a]] -> [([a], [a])]
    separateTrainAndValidLists = fmap (\(t,v)->(t [], v)) . separate' ([]++)

    separate' :: ([a]->[a]) -> [[a]] -> [([a]->[a], [a])]
    separate' acc [x]    = [(acc, x)]
    separate' acc (x:xs) = (acc . (concat xs++), x) : separate' (acc . (x++)) xs
    separate' _   []     = error "this should never happen, the size of the partition is too big to even hold an element"

saveScores :: FilePath -> [NetScore] -> IO ()
saveScores logsPath scores = writeFile logsPath (headScores<>"\n"<>scoresStr) -- TODO: catch IO Errors!
  where
    scoresStr :: String
    scoresStr = mconcat . fmap netScore2String $ zip [(0::Int)..] scores

    netScore2String (i, NetScore{..}) = show i <> "\t"
                                     <> show trainClassError <> "\t"
                                     <> show trainingError   <> "\t"
                                     <> show testClassError  <> "\t"
                                     <> show gradientNorm    <> "\n"

    headScores = "epoch\ttrain_classification_error\ttrain_error\ttest_classification_error\tgradient_norm"

-- Reading data and running code

filenameDataset :: String
filenameDataset = "./msd_genre_dataset.txt"

getSongs :: String -> IO [Song]
getSongs filename = do
  --putStrLn "Reading data file ..." -- <- this doesn't work by the lazy nature of haskell
  file <- Text.readFile filename
  --let listLines = drop 10 $ Data.Text.lines file
  let listLines = Data.Text.lines file
      all_songs = mapMaybe line2song listLines
  return $ filter inRightGenre all_songs

 where inRightGenre song = let g = view genre song
                            in g == "dance and electronica" || g == "jazz and blues"

labelSong :: String -> Double
labelSong "dance and electronica" = 1
labelSong _                       = 0

data StoppingCondition = StoppingCondition { -- stop training if either:
  maxEpochs           :: Int,           -- maximum total epochs has been reached
  gradientSmallerthan :: (Double, Int), -- gradient is smaller than `d` for at least `i` consecutive epochs
  stuckInARange       :: (Double, Int)  -- or, the error hasn't changed more than `d` in the last `i` consecutive epochs
  }

data TypeOfModelToTrain = LogitRegModel | OneHiddenLayerModel Integer | ArbitraryNNModel Int

data ModelsParameters =
  ModelsParameters
    Float -- Size of test set
    Int   -- Size of batch
    Bool  -- Normalize data
    StoppingCondition
    LearningParameters -- Parameters for Gradient Descent
    (Maybe FilePath) -- Load path
    (Maybe FilePath) -- Save path
    (Maybe FilePath) -- Logs path
    (Maybe Double)   -- K-fold size of validation set
    TypeOfModelToTrain

modelsParameters :: Parser ModelsParameters
modelsParameters =
  ModelsParameters <$> option auto (long "test-set"   <> short 't' <> value 0.15)
                   <*> option auto (long "batch-size" <> short 'b' <> value 40)
                   <*> option auto (long "normalize"  <> value True)
                   <*> ( StoppingCondition
                         <$> option auto (long "epochs" <> short 'e' <> value 300)
                         <*> option auto (long "gradient-smaller" <> short 'g' <> value (0.085, 3))
                         <*> option auto (long "max-error-change" <> value (0.004, 15))
                       )
                   <*> (LearningParameters
                       <$> option auto (long "train_rate" <> short 'r' <> value 0.01)
                       <*> option auto (long "momentum" <> value 0.9)
                       <*> option auto (long "l2" <> value 0) -- 0.0005
                       )
                   <*> optional (strOption (long "load"))
                   <*> optional (strOption (long "save"))
                   <*> optional (strOption (long "logs"))
                   <*> optional (option auto (long "k-fold-val-size"))
                   <*> (   flag' LogitRegModel (long "logit-reg")
                       <|> (OneHiddenLayerModel <$> option auto (long "one-hidden-layer" <> short 'h'))
                       <|> (ArbitraryNNModel    <$> option auto (long "arbitrary-nn"))
                       )

addLog :: KnownNat n => R n -> R n
addLog a = signum a * log (abs a+1)

discreteToOneHotVector :: R 3 -> R 21
discreteToOneHotVector featDiscre = (fromMaybe (konst 0) $ hotvector (round a)::R 8)
                                  # (fromMaybe (konst 0) $ hotvector (round b)::R 12)
                                  & c
  where (a,b,c) = case toList $ unwrap featDiscre of
                    [a',b',c'] -> (a',b',c')
                    _          -> error $  "This is weird, there should be only three discrete "
                                        <> "features (R 3), aka. this never happens but the compiler "
                                        <> "asks me to give an exahustive list for pattern matching, doh'"

takeWhileCond :: ([a] -> [b]) -> (b -> Bool) -> [a] -> [a]
takeWhileCond toCond stopCond xs = fmap snd . takeWhile (stopCond . fst) $ zip (toCond xs) xs

data NetScore = NetScore {
  trainClassError :: Double, -- Classification training error
  trainingError   :: Double, -- (optimization) training error
  testClassError  :: Double, -- Classification testing error
  gradientNorm    :: Double  -- accumulated norm of the gradients on an epoch
} deriving (Show)

netScore :: (Last shapes ~ 'D1 1, Foldable t)
         => t (S (Head shapes), S ('D1 1))
         -> t (S (Head shapes), S ('D1 1))
         -> (Double, Network layers shapes)
         -> NetScore
netScore trainSet testSet (gradNorm, nn) =
  let (trainCE, trainE) = binaryNetError nn trainSet
      (testCE, _)       = binaryNetError nn testSet
   in
      NetScore { trainClassError = trainCE,
                 trainingError   = trainE,
                 testClassError  = testCE,
                 gradientNorm    = gradNorm }

takeWhileCondFunc :: Show a => StoppingCondition -> [(a, NetScore)] -> [(a, NetScore)]
takeWhileCondFunc (StoppingCondition epochs (grad, giters) (stuck, siters)) =
 -- trace ("hi" <> show (head (takeWhileCond maxErrorChangeInSIters (<stuck) xs))) .
  take (epochs+1) . takeWhileCond consecutiveSmallerThanGrad (<giters)
                  . takeWhileCond maxErrorChangeInSIters (>stuck)
  where
    consecutiveSmallerThanGrad :: [(a, NetScore)] -> [Int]
    consecutiveSmallerThanGrad = (0:) . countConsecutive (0::Int) . fmap (trainClassError . snd)
      where
        countConsecutive _    []     = []
        countConsecutive cond (x:xs) =
          let condVal = if x<grad then cond+1 else 0
           in condVal : countConsecutive condVal xs

    maxErrorChangeInSIters :: [(a, NetScore)] -> [Double]
    maxErrorChangeInSIters = (replicate siters (stuck+1) <>) . fmap findMaxChange . groupOnSIters . fmap (trainClassError . snd)
      where groupOnSIters :: [a] -> [[a]]
            groupOnSIters []         = []
            groupOnSIters xs@(_:xs') = take siters xs : groupOnSIters xs'
            findMaxChange :: [Double] -> Double
            findMaxChange xs = maximum xs - minimum xs

loadModel (Just loadFile) _       = netLoad loadFile
loadModel Nothing         seedNet = return $ evalRand randomNetwork seedNet

main :: IO ()
main = do
  mparams@(ModelsParameters _ _ _ _ _ load _ _ _ modelType) <- execParser (info (modelsParameters <**> helper) idm)

  let (seedNet, randSeed) = SR.split (mkStdGen 487239842)

  case modelType of
    LogitRegModel ->
      (loadModel load seedNet :: IO LogisticRegression) >>= mainWithRandmNet mparams randSeed

    OneHiddenLayerModel sizeHidden ->
      let (singSizeHidden :: SomeSing Nat) = toSing sizeHidden
      in case singSizeHidden of
           SomeSing (SNat :: Sing n) ->
             (loadModel load seedNet :: IO (OneHiddenLayer n)) >>=
             mainWithRandmNet mparams randSeed

    ArbitraryNNModel arbModelNum ->
      case arbModelNum of
        1 -> (loadModel load seedNet :: IO DiscreteSepContinuousNet) >>= mainWithRandmNet mparams randSeed
        2 -> (loadModel load seedNet :: IO HugeNetwork) >>= mainWithRandmNet mparams randSeed
        _ -> putStrLn $ "Sorry but there is no arbitrary neural net number " <> show arbModelNum

mainWithRandmNet :: (Head shapes ~ 'D1 75,
                     Last shapes ~ 'D1 1,
                     Show (Network layers shapes),
                     Serialize (Network layers shapes),
                     GradNorm (Gradients layers),
                     Num (Gradients layers),
                     RandomGen g)
                  => ModelsParameters -> g -> Network layers shapes -> IO ()
mainWithRandmNet mparams randSeed net0 = do
  let (seedShuffle, seedTraining) = SR.split randSeed
      ModelsParameters testPerc batchSize norm stopCond rate _ save logs kfoldValPer modelType = mparams
      modelName = case modelType of
                    LogitRegModel                  -> fromMaybe "logitReg"  logs
                    OneHiddenLayerModel sizeHidden -> fromMaybe "hidden"    logs <> "-" <> show sizeHidden
                    ArbitraryNNModel arbModelNum   ->
                      fromMaybe "arbitrary" logs <> "-" <> case arbModelNum of
                                                             1 -> "separated_discrete_continuous"
                                                             2 -> "huge_network"
                                                             _ -> error $ "There is no arbitrary Neural Net model #" <> show arbModelNum

  putStrLn $ "Model name: " <> modelName

  let normFun :: KnownNat m => [R m] -> [R m]
      normFun = if norm then normalize else id

  -- Loading songs
  songsRaw <- fmap (song2TupleRn labelSong) <$> getSongs filenameDataset
  let songsDiscrete :: [R 3]
      songsDiscrete = fmap (fst . fst) songsRaw
      songsDisOneHotVector :: [R 21]
      songsDisOneHotVector = fmap discreteToOneHotVector songsDiscrete
      --songsFloat    :: [R 27]
      songsFloat    = fmap (snd . fst) songsRaw
      --songsLog      :: [R 27]
      songsLog      = fmap addLog songsFloat
      --songsLabel    :: [R 1]
      songsLabel    = fmap snd songsRaw
      --songs :: [(S ('D1 57), S ('D1 1))]
      songs = (\(l1,l2,l3,o)->(S1D (l1#l2#l3), S1D o)) <$> zip4 songsDisOneHotVector (normFun songsFloat) (normFun songsLog) songsLabel

  --print . toList . unwrap $ foldl1' (zipWithVector max) songsDiscrete
  --print . toList . unwrap $ foldl1' (zipWithVector min) songsDiscrete

  let sizeSongs = length songsRaw
      -- Shuffling songs
      shuffledSongs       = evalRand (shuffle songs) seedShuffle
      (testSet, trainSet) = splitAt (round $ fromIntegral sizeSongs * testPerc) shuffledSongs

  {-
   -case net0 of
   -  x@(FullyConnected (FullyConnected' i o) _) :~> _ -> do
   -    let subnet = x :~> NNil :: Network '[FullyConnected 75 n] '[ 'D1 75, 'D1 n ]
   -    print $ runNet subnet (S1D 0) -- with this we get the biases
   -    --print . sqrt . normSquared $ net0
   -    print i -- biases for each output neuron
   -    print o -- weights between neurons
   -}

  -- Pretraining scores
  --print $ head shuffledSongs
  putStrLn $ "Total Size: " <> show sizeSongs
  putStrLn $ "Test Size:  " <> show (length testSet)
  putStrLn $ "Batch Size: " <> show batchSize

  case kfoldValPer of
    Nothing -> trainNet net0 trainSet testSet rate batchSize seedTraining stopCond logs save
    Just valPer -> do
      let valSize             = floor $ fromIntegral sizeSongs * valPer
          (parts, kfoldSets)  = createKFoldPartitions valSize trainSet
      putStrLn $ "Validation Size: " <> show valSize
      putStrLn $ "Total parts for kfold: " <> show parts

      forM_ (zip [(1::Int)..] kfoldSets) $ \(kfoldNum, (trainSet', valSet)) -> do
        let dirSaveNet = "kfold-val" <> show valPer <> "/" <> modelName <> "/" <> "kfoldnet" <> "-" <> show kfoldNum
            logsKfoldN = dirSaveNet <> ".txt"
            saveKfoldN = dirSaveNet <> "/" <> fromMaybe "nnet-hidden" save
        putStrLn ""
        createDirectoryIfMissing True dirSaveNet -- TODO: catch exceptions
        putStrLn $ "Results to be saved in " <> logsKfoldN
        putStrLn $ "Nets to be saved in " <> saveKfoldN
        trainNet net0 trainSet' valSet rate batchSize seedTraining stopCond (Just logsKfoldN) (Just saveKfoldN)

trainNet :: (Last shapes ~ 'D1 1,
             RandomGen g,
             Num (Gradients layers),
             GradNorm (Gradients layers),
             Show (Network layers shapes),
             Serialize (Network layers shapes))
         => Network layers shapes
         -> [(S (Head shapes), S (Last shapes))]
         -> [(S (Head shapes), S ('D1 1))]
         -> LearningParameters
         -> Int
         -> g
         -> StoppingCondition
         -> Maybe FilePath
         -> Maybe FilePath
         -> IO ()
trainNet net0 trainSet testSet rate batchSize seedTraining stopCond logs save = do
  -- Training Net
  let netsInf = (read "Infinity", net0) : evalRand (trainOnBatchesEpochs net0 rate trainSet batchSize) seedTraining
      netsScoresInf = fmap (netScore trainSet testSet) netsInf

      (nets, netsScores) = unzip . takeWhileCondFunc stopCond $ zip (fmap snd netsInf) netsScoresInf
      --net = last nets

  putStrLn $ "Epoch"
           <> "\tTraining classification error"
           <> "\tTraining error"
           <> "\tTesting error"
           <> "\tGradNorm"

  -- Showing results of training net
  forM_ (zip3 [(0::Integer)..] netsScores nets) $ \(epoch, NetScore trainCE trainE testE gradNorm, currentNet) -> do
    putStrLn $ show epoch
             <> "\t" <> show trainCE
             <> "\t" <> show trainE
             <> "\t" <> show testE
             <> "\t" <> show gradNorm
    -- saving current trained net
    case save of
      Just saveFile -> let saveFileBin = saveFile <> "-e_" <> show epoch
                                                  <> "-trainCE_" <> show trainCE
                                                  <> "-testE_"   <> show testE
                                                  <> ".bin"
                        in B.writeFile saveFileBin $ runPut (put currentNet)
      Nothing -> return ()

  case logs of
    Just logsPath -> saveScores logsPath netsScores
    Nothing       -> return ()

  --case save of
  --  Just saveFile -> B.writeFile saveFile $ runPut (put net)
  --  Nothing       -> return ()

max_stars_repo_name: helq/haskell-binary-classification

max_issues_repo_name: helq/haskell-binary-classification

max_forks_repo_name: helq/haskell-binary-classification


text:
{-# LANGUAGE BangPatterns #-}

module BBrotRender(showCells, render) where

import Codec.Picture
import Control.Monad
import Data.Aeson(decode)
import Data.Array.IO
import Data.Array.Unboxed
import qualified Data.ByteString.Lazy as BS
import Data.Complex
import Data.Maybe
import Data.List
import Data.Word(Word32, Word8)
import System.Console.CmdArgs(whenNormal)
import System.Exit
import System.IO
import Text.Read(readMaybe)

import BBrotCompute
import BBrotConf
import BBrotSelection

xmin, xmax, ymin, ymax :: Double
xmin = realPart loCorner
xmax = realPart hiCorner
ymin = imagPart loCorner
ymax = imagPart hiCorner

xrange, yrange :: Double
xrange = xmax - xmin
yrange = ymax - ymin

orbs :: Complex Double -> [Complex Double]
orbs !(x0 :+ y0) = unfoldr f (x0, y0)
  where f :: (Double, Double) -> Maybe (Complex Double, (Double, Double))
        f (!x, !y) =
          let !x2 = x * x
              !y2 = y * y
              !newx = (x2 - y2 + x0)
              !newy = (2 * x * y + y0)
          in
           if x2 + y2 > 4
           then Nothing
           else Just (newx :+ newy, (newx, newy))

inWindow :: Complex Double -> Bool
inWindow z = x >= xmin && x < xmax && y >= ymin && y < ymax
  where x = realPart z
        y = imagPart z

norm :: (Double -> Double) -> Word32 -> Word32 -> Word32 -> Word8
norm curveFunc minV maxV cnt = v
  where t = fromIntegral (cnt - minV) / fromIntegral (maxV - minV)
        u = curveFunc $ min 1.0 (2.0 * t)
        v = floor $ u * 255 :: Word8

gray :: Word8 -> PixelRGB8
gray x = PixelRGB8 x x x

reddish :: Word8 -> PixelRGB8
reddish x = PixelRGB8 r g 0
  where xr = min x 127
        xg = min (x - xr) 127
        r = 2 * xr
        g = 2 * xg

flames :: Word8 -> PixelRGB8
flames x = PixelRGB8 r g b
  where xr = min 120 x
        xg = min 120 (x - xr)
        xb = min 15 (x - xr - xg)
        r = xr * 2 + xr `div` 10
        g = xg * 2 + xg `div` 10
        b = xb * 17

toPixel :: (Double -> Double) -> Word32 -> Word32 -> (Word8 -> PixelRGB8) -> Word32 -> PixelRGB8
toPixel curveFunc smallest biggest scheme = scheme . norm curveFunc smallest biggest

plotPix :: IOUArray (Int, Int) Word32 -> (Int, Int) -> IO ()
plotPix !img (!x, !y) = do
  v <- readArray img (x, y)
  writeArray img (x, y) (v + 1)

rel :: Double -> Double -> Double -> Double
rel !lo !realrange !v = (v - lo) / realrange

toPix :: Int -> Double -> Double -> Double -> Int
toPix !pixrange !realmin !realrange !a =
  floor $ (fromIntegral pixrange) * (rel realmin realrange a)

toImgCoords :: Int -> Int -> Complex Double -> (Int, Int)
toImgCoords !xres !yres !z =
    (toPix xres xmin xrange $ realPart z,
     toPix yres ymin yrange $ imagPart z)

toPlaneCoords :: Int -> Int -> Int -> Int -> Complex Double
toPlaneCoords !xres !yres !i !j = x :+ y
    where x = xmin + xrange * fromIntegral i / fromIntegral xres
          y = ymin + yrange * fromIntegral j / fromIntegral yres

emptyPS :: PointSelection
emptyPS = PointSelection { pointList   = []
                         , commandLine = Nothing
                         , randGen     = Nothing
                         , timeStamp   = Nothing
                         }

loadPointsJson :: String -> IO (Maybe PointSelection)
loadPointsJson filepath = do
  contents <- BS.readFile filepath
  return (decode contents :: Maybe PointSelection)

loadPointsComplex :: String -> IO (Maybe PointSelection)
loadPointsComplex filepath = do
  contents <- readFile filepath
  let zs = concatMap g $ lines contents
        where g = maybeToList . readMaybe :: String -> [Complex Double]
      ps = map f zs
        where f (x :+ y) = BBPoint x y 1
  return $ if null ps
           then Nothing
           else Just (emptyPS { pointList = ps })

render :: BBrotConf -> IO ()
render conf = do
  -- Load points of interest, render their orbits into a PNG image.
  whenNormal $ putStrLn $ "Loading cache " ++ icachepath conf ++ " ..."
  maybePS <- if isComplex conf
             then loadPointsComplex $ icachepath conf
             else loadPointsJson $ icachepath conf
  when (isNothing maybePS) $ do
    hPutStrLn stderr $ "failed to parse " ++ icachepath conf
    unless (isComplex conf) $
      hPutStrLn stderr "try with flag -z"
    exitFailure

  let psel = pointList $ fromMaybe emptyPS maybePS
      selected = map (\(BBPoint x y _) -> x :+ y) psel
      orbits = concatMap orbs selected
      result = filter inWindow $ if dontRender conf
                                 then selected
                                 else orbits
  whenNormal $ putStrLn $ "selected points: " ++ show (length selected)

  let xres = xpixels conf
      yres = ypixels conf
      coords = map (toImgCoords xres yres) result
  img <- newArray ((0, 0), (xres - 1, yres - 1)) 0 :: IO (IOUArray (Int, Int) Word32)
  mapM_ (plotPix img) coords
  whenNormal $ putStrLn "done plotting"

  img2 <- freeze img :: IO (UArray (Int, Int) Word32)
  let values = elems img2
      !v0 = head values
      (!total, !smallest, !biggest) = foldl' f (0, v0, v0) values
        where f (!x, !y, !z) !a = (a + x, min a y, max a z)

  whenNormal $ do
    putStrLn $ "img points: " ++ show total
    putStrLn $ "value range: " ++ show smallest ++ "-" ++ show biggest

  let outfile = fromMaybe defPath (imagepath conf)
        where defPath = icachepath conf ++ ".png"
  whenNormal $ putStrLn $ "Writing " ++ outfile ++ " ..."
  let colorScheme = case palette conf of
        Flames -> flames
        Gray -> gray
        Reddish -> reddish
      curveFunc = case curve conf of
        Line -> id
        Root -> sqrt
        Square -> (**2)
      pixFunc = toPixel curveFunc smallest biggest colorScheme
      renderer i j = pixFunc $ img2!(i,j)
  writePng outfile $ generateImage (flip $ renderer) yres xres

showCells :: BBrotConf -> IO ()
showCells conf = do
  -- Compute selection of cells (squares in the complex plane) that
  -- are close to the edge of the Mandelbrot set, and render them as
  -- animated GIF.
  let step = gridStep conf
      bailout = maxIters conf
      xres = 1000
      yres = 1000
      !cells = selectCells step bailout
  whenNormal $ putStrLn $ "cell count: " ++ show (length cells)

  let red   = PixelRGB8 255 0 0
      black = PixelRGB8 0 0 0
      grey  = PixelRGB8 64 64 64

  whenNormal $ putStrLn "rendering cells"
  cellMap <- newArray ((0, 0), (xres - 1, yres - 1)) False :: IO (IOUArray (Int, Int) Bool)
  forM_ cells $ \(x, y) -> do
    let cornerA = (x - step / 2) :+ (y - step / 2)
        cornerB = (x + step / 2) :+ (y + step / 2)
        (imin, jmin) = toImgCoords xres yres cornerA
        (imax, jmax) = toImgCoords xres yres cornerB
        coords = [ (i, j) | i <- [imin..imax], j <- [jmin..jmax],
                                 i >= 0 && i < xres,
                                 j >= 0 && j < yres ]
    forM_ coords $ \(i, j) ->
      writeArray cellMap (i, j) True

  whenNormal $ putStrLn "rendering mandel"
  let imgMandel = generateImage (flip $ mandelRenderer) yres xres
      inMandelbrotSet z = inSet 0 bailout z
      mandelRenderer i j = if inMandelbrotSet $ toPlaneCoords xres yres i j
                           then grey else black

  cellMapPure <- freeze cellMap :: IO (UArray (Int, Int) Bool)
  let imgCells = generateImage (flip $ cellRenderer) yres xres
      cellRenderer i j = if cellMapPure!(i,j)
                         then red else mandelRenderer i j

  whenNormal $ putStrLn $ "writing " ++ animpath conf ++ " ..."
  case writeGifAnimation (animpath conf) 100 LoopingForever [imgMandel, imgCells] of
    Left err -> putStrLn $ "error generating gif animation: " ++ err
    Right action -> action

ext: hs

max_stars_repo_name: saffroy/buddhabrot

max_issues_repo_name: saffroy/buddhabrot

max_forks_repo_name: saffroy/buddhabrot


text:
{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TemplateHaskell     #-}

module Test.Grenade.Loss where

import           Hedgehog
import qualified Hedgehog.Gen                 as Gen
import qualified Hedgehog.Range               as Range

import qualified Numeric.LinearAlgebra.Static as NLA

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

genShapes :: Gen [S ('D1 10)]
genShapes =
  let len = Range.linear 0 10
  in  Gen.list len (S1D <$> randomVectorNormalised)

prop_quadratic :: Property
prop_quadratic = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry quadratic) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> 0.5 * (sum $ zipWith (\a b -> (a - b) ^ (2 :: Int)) x y)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_quadratic' :: Property
prop_quadratic' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = quadratic'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (-)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_crossEntropy :: Property
prop_crossEntropy = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry crossEntropy) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> -1 * (sum $ zipWith (\a b -> b * (log a) + (1 - b) * (log (1 - a))) x y)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_crossEntropy' :: Property
prop_crossEntropy' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = crossEntropy'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = (zipWith (\a b -> (a - b) / ((1 - a) * a)))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_exponential :: Property
prop_exponential = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let t  = 0.1
  let costs  = map (uncurry (exponential t)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> t * (exp (1/t * sum (zipWith (\a b -> (a - b) ^ (2 :: Int)) x y)))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_exponential' :: Property
prop_exponential' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let t = 0.1
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = exponential' t
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let g      = \x y -> map ((t * (exp (1/t * sum (zipWith (\a b -> (a - b) ^ (2 :: Int)) x y)))) *)
  let f      = \x y -> map ((2 / t) *) $ zipWith (-) x y
  let costs' = zipWith (\x y -> g x y (f x y)) xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_hellinger :: Property
prop_hellinger = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry hellinger) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> (1 / (sqrt 2)) * (sum $ zipWith (\a b -> (sqrt a - sqrt b) ^ (2 :: Int)) x y)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_hellinger' :: Property
prop_hellinger' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = hellinger'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (\a b -> ((sqrt a) - (sqrt b)) / ((sqrt 2) * (sqrt a)))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_kullbackLeibler :: Property
prop_kullbackLeibler = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry kullbackLeibler) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> sum $ zipWith (\a b -> b * log (b / a)) x y
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_kullbackLeibler' :: Property
prop_kullbackLeibler' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = kullbackLeibler'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (\a b -> -(b / a))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_genKullbackLeibler :: Property
prop_genKullbackLeibler = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry genKullbackLeibler) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> (sum $ zipWith (\a b -> b * log (b / a)) x y) - (sum y) + (sum x)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_genKullbackLeibler' :: Property
prop_genKullbackLeibler' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = genKullbackLeibler'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (\a b -> ((a - b) / a))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_itakuraSaito :: Property
prop_itakuraSaito = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry itakuraSaito) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> sum $ zipWith (\a b -> (b/a) - (log (b/a)) - 1) x y
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_itakuraSaito' :: Property
prop_itakuraSaito' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = itakuraSaito'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (\a b -> ((a - b) / (a * a)))
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

prop_categoricalCrossEntropy :: Property
prop_categoricalCrossEntropy = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let costs  = map (uncurry categoricalCrossEntropy) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = \x y -> -1 * (sum $ zipWith (\a b -> b * log a) x y)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListTo` costs'

prop_categoricalCrossEntropy' :: Property
prop_categoricalCrossEntropy' = property $ do
  xs <- forAll genShapes
  ys <- forAll genShapes
  let zs = zip xs ys
  let (LossFunction (l :: (S ('D1 10) -> S ('D1 10) -> S ('D1 10)))) = categoricalCrossEntropy'
  let costs  = map (extractVec . (uncurry l)) zs
  let xs'    = map extractVec xs
  let ys'    = map extractVec ys
  let f      = zipWith (\a b -> -b / a)
  let costs' = zipWith f xs' ys'
  costs `isSimilarListOfListsTo` costs'

exampleShape :: (S ('D1 3), S ('D1 3))
exampleShape = (gen predicted, gen true)
    where
      gen       = S1D . NLA.vector
      predicted = [0.1, 0.8, 0.1]
      true      = [0, 0, 1]

prop_categoricalCrossEntropy_matchesActualValues :: Property
prop_categoricalCrossEntropy_matchesActualValues = property $ do
    let calculated = (uncurry categoricalCrossEntropy) exampleShape :: RealNum
    diff (abs (calculated - expectedResult)) (<) 0.01
    where
      expectedResult = 2.303

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

max_stars_repo_name: th-char/grenade

max_issues_repo_name: th-char/grenade

max_forks_repo_name: th-char/grenade


text:
-- https://github.com/minoki/my-atcoder-solutions
{-# LANGUAGE BangPatterns          #-}
{-# LANGUAGE DataKinds             #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE NoStarIsType          #-}
{-# LANGUAGE ScopedTypeVariables   #-}
{-# LANGUAGE TypeApplications      #-}
{-# LANGUAGE TypeFamilies          #-}
{-# LANGUAGE TypeOperators         #-}
import           Control.Exception           (assert)
import           Control.Monad
import           Data.Bits
import qualified Data.ByteString.Char8       as BS
import           Data.Char                   (isSpace)
import           Data.Coerce
import           Data.Complex
import           Data.Int                    (Int64)
import           Data.List                   (unfoldr)
import qualified Data.Vector.Generic         as G
import qualified Data.Vector.Unboxing         as U
import qualified Data.Vector.Unboxing.Mutable as UM
import           GHC.TypeNats                (type (+), type (*), KnownNat, Nat, type (^),
                                              natVal)

type R1 = IntMod (5 * 2^25 + 1)
type R2 = IntMod (7 * 2^26 + 1)
type R3 = IntMod (3 * 2^18 + 1)
type R4 = IntMod (7 * 2^20 + 1)
type R5 = IntMod (483 * 2^21 + 1)

order :: (Num a, Eq a) => a -> Int
order !x = go 1 x
  where
    go !n 1 = n
    go !n y = go (n + 1) (x * y)

order' :: (Num a, Eq a) => Int -> a -> Int
order' !m !x = go 1 x
  where
    go !n 1 = n
    go !n y | n > m = m + 1
    go !n y = go (n + 1) (x * y)

findPrimitiveNthRoot :: (Num a, Eq a) => Int -> a
findPrimitiveNthRoot n = head [ x | k <- [1..], let x = fromInteger k, order' n x == n ]

main = do
  -- print (findPrimitiveNthRoot 2)
  -- print (findPrimitiveNthRoot (2^10))
  print (findPrimitiveNthRoot (2^25) :: R1) -- 17 mod 5 * 2^25 + 1
  print (findPrimitiveNthRoot (2^26) :: R2) -- 30 mod 7 * 2^26 + 1
  print (findPrimitiveNthRoot (2^18) :: R3) -- 5 mod 3 * 2^18 + 1
  print (findPrimitiveNthRoot (2^20) :: R4) -- 5 mod 7 * 2^20 + 1
  print (findPrimitiveNthRoot (2^21) :: R5) -- 198 mod 483 * 2^21 + 1

--
-- Fast Fourier Transform (FFT)
--

halve :: G.Vector vec a => vec a -> vec a
halve v = let n = G.length v
          in G.generate (n `quot` 2) $ \j -> v G.! (j * 2)

fft :: forall vec a. (Num a, G.Vector vec a)
    => [vec a] -- ^ For a primitive n-th root of unity @u@, @iterate halve [1,u,u^2 .. u^(n-1)]@
    -> vec a -- ^ a polynomial of length n (= 2^k for some k)
    -> vec a
fft (u:u2) f | n == 1 = f
             | otherwise = let !n2 = n `quot` 2
                               r0, r1', t0, t1' :: vec a
                               r0 = G.generate n2 $ \j -> (f G.! j) + (f G.! (j + n2))
                               r1' = G.generate n2 $ \j -> ((f G.! j) - (f G.! (j + n2))) * u G.! j
                               !t0 = fft u2 r0
                               !t1' = fft u2 r1'
                           in G.generate n $ \j -> if even j then t0 G.! (j `quot` 2) else t1' G.! (j `quot` 2)
  where n = G.length f

mulFFT :: U.Vector Int -> U.Vector Int -> U.Vector Int
mulFFT !f !g = let n' = U.length f + U.length g - 2
                   k = finiteBitSize n' - countLeadingZeros n'
                   !_ = assert (n' < 2^k) ()
                   n = bit k
                   u :: U.Vector (Complex Double)
                   u = U.generate n $ \j -> cis (fromIntegral j * (2 * pi / fromIntegral n))
                   us = iterate halve u
                   f' = U.generate n $ \j -> if j < U.length f then
                                               fromIntegral (f U.! j)
                                             else
                                               0
                   g' = U.generate n $ \j -> if j < U.length g then
                                               fromIntegral (g U.! j)
                                             else
                                               0
                   f'' = fft us f'
                   g'' = fft us g'
                   fg = U.generate n $ \j -> (f'' U.! j) * (g'' U.! j)
                   fg' = fft (map (U.map conjugate) us) fg
               in U.generate n $ \j -> round (realPart (fg' U.! j) / fromIntegral n)

--
-- Univariate polynomial
--

newtype Poly vec a = Poly { coeffAsc :: vec a } deriving Eq

normalizePoly :: (Eq a, Num a, G.Vector vec a) => vec a -> vec a
normalizePoly v | G.null v || G.last v /= 0 = v
                | otherwise = normalizePoly (G.init v)

addPoly :: (Eq a, Num a, G.Vector vec a) => vec a -> vec a -> vec a
addPoly v w = case compare n m of
                LT -> G.generate m $ \i -> if i < n
                                           then v G.! i + w G.! i
                                           else w G.! i
                GT -> G.generate n $ \i -> if i < m
                                           then v G.! i + w G.! i
                                           else v G.! i
                EQ -> normalizePoly $ G.zipWith (+) v w
  where n = G.length v
        m = G.length w

subPoly :: (Eq a, Num a, G.Vector vec a) => vec a -> vec a -> vec a
subPoly v w = case compare n m of
                LT -> G.generate m $ \i -> if i < n
                                           then v G.! i - w G.! i
                                           else negate (w G.! i)
                GT -> G.generate n $ \i -> if i < m
                                           then v G.! i - w G.! i
                                           else v G.! i
                EQ -> normalizePoly $ G.zipWith (-) v w
  where n = G.length v
        m = G.length w

naiveMulPoly :: (Num a, G.Vector vec a) => vec a -> vec a -> vec a
naiveMulPoly v w = G.generate (n + m - 1) $
                   \i -> sum [(v G.! (i-j)) * (w G.! j) | j <- [max (i-n+1) 0..min i (m-1)]]
  where n = G.length v
        m = G.length w

doMulP :: (Eq a, Num a, G.Vector vec a) => Int -> vec a -> vec a -> vec a
doMulP n !v !w | n <= 16 = naiveMulPoly v w
doMulP n !v !w
  | G.null v = v
  | G.null w = w
  | G.length v < n2 = let (w0, w1) = G.splitAt n2 w
                          u0 = doMulP n2 v w0
                          u1 = doMulP n2 v w1
                      in G.generate (G.length v + G.length w - 1)
                         $ \i -> case () of
                                   _ | i < n2     ->  u0 `at` i
                                     | i < n      -> (u0 `at` i) + (u1 `at` (i - n2))
                                     | i < n + n2 -> (u1 `at` (i - n2))
  | G.length w < n2 = let (v0, v1) = G.splitAt n2 v
                          u0 = doMulP n2 v0 w
                          u1 = doMulP n2 v1 w
                      in G.generate (G.length v + G.length w - 1)
                         $ \i -> case () of
                                   _ | i < n2     ->  u0 `at` i
                                     | i < n      -> (u0 `at` i) + (u1 `at` (i - n2))
                                     | i < n + n2 -> (u1 `at` (i - n2))
  | otherwise = let (v0, v1) = G.splitAt n2 v
                    (w0, w1) = G.splitAt n2 w
                    v0_1 = v0 `addPoly` v1
                    w0_1 = w0 `addPoly` w1
                    p = doMulP n2 v0_1 w0_1
                    q = doMulP n2 v0 w0
                    r = doMulP n2 v1 w1
                    -- s = (p `subPoly` q) `subPoly` r -- p - q - r
                    -- q + s*X^n2 + r*X^n
                in G.generate (G.length v + G.length w - 1)
                   $ \i -> case () of
                             _ | i < n2     ->   q `at` i
                               | i < n      -> ((q `at` i) + (p `at` (i - n2))) - ((q `at` (i - n2)) + (r `at` (i - n2)))
                               | i < n + n2 -> ((r `at` (i - n)) + (p `at` (i - n2))) - ((q `at` (i - n2)) + (r `at` (i - n2)))
                               | otherwise  ->   r `at` (i - n)
  where n2 = n `quot` 2
        at :: (Num a, G.Vector vec a) => vec a -> Int -> a
        at v i = if i < G.length v then v G.! i else 0
{-# INLINE doMulP #-}

mulPoly :: (Eq a, Num a, G.Vector vec a) => vec a -> vec a -> vec a
mulPoly !v !w = let k = ceiling ((log (fromIntegral (max n m)) :: Double) / log 2) :: Int
                in doMulP (2^k) v w
  where n = G.length v
        m = G.length w
{-# INLINE mulPoly #-}

zeroPoly :: (G.Vector vec a) => Poly vec a
zeroPoly = Poly G.empty

constPoly :: (Eq a, Num a, G.Vector vec a) => a -> Poly vec a
constPoly 0 = Poly G.empty
constPoly x = Poly (G.singleton x)

scalePoly :: (Eq a, Num a, G.Vector vec a) => a -> Poly vec a -> Poly vec a
scalePoly a (Poly xs)
  | a == 0 = zeroPoly
  | otherwise = Poly $ G.map (* a) xs

valueAtPoly :: (Num a, G.Vector vec a) => Poly vec a -> a -> a
valueAtPoly (Poly xs) t = G.foldr' (\a b -> a + t * b) 0 xs

instance (Eq a, Num a, G.Vector vec a) => Num (Poly vec a) where
  (+) = coerce (addPoly :: vec a -> vec a -> vec a)
  (-) = coerce (subPoly :: vec a -> vec a -> vec a)
  negate (Poly v) = Poly (G.map negate v)
  (*) = coerce (mulPoly :: vec a -> vec a -> vec a)
  fromInteger = constPoly . fromInteger
  abs = undefined; signum = undefined

divModPoly :: (Eq a, Fractional a, G.Vector vec a) => Poly vec a -> Poly vec a -> (Poly vec a, Poly vec a)
divModPoly f g@(Poly w)
  | G.null w = error "divModPoly: divide by zero"
  | degree f < degree g = (zeroPoly, f)
  | otherwise = loop zeroPoly (scalePoly (recip b) f)
  where
    g' = toMonic g
    b = leadingCoefficient g
    -- invariant: f == q * g + scalePoly b p
    loop q p | degree p < degree g = (q, scalePoly b p)
             | otherwise = let q' = Poly (G.drop (degree' g) (coeffAsc p))
                           in loop (q + q') (p - q' * g')

    toMonic :: (Fractional a, G.Vector vec a) => Poly vec a -> Poly vec a
    toMonic f@(Poly xs)
      | G.null xs = zeroPoly
      | otherwise = Poly $ G.map (* recip (leadingCoefficient f)) xs

    leadingCoefficient :: (Num a, G.Vector vec a) => Poly vec a -> a
    leadingCoefficient (Poly xs)
      | G.null xs = 0
      | otherwise = G.last xs

    degree :: G.Vector vec a => Poly vec a -> Maybe Int
    degree (Poly xs) = case G.length xs - 1 of
      -1 -> Nothing
      n  -> Just n

    degree' :: G.Vector vec a => Poly vec a -> Int
    degree' (Poly xs) = case G.length xs of
      0 -> error "degree': zero polynomial"
      n -> n - 1

-- 組立除法
-- second constPoly (divModByDeg1 f t) = divMod f (Poly (G.fromList [-t, 1]))
divModByDeg1 :: (Eq a, Num a, G.Vector vec a) => Poly vec a -> a -> (Poly vec a, a)
divModByDeg1 f t = let w = G.postscanr (\a b -> a + b * t) 0 $ coeffAsc f
                   in (Poly (G.tail w), G.head w)

--
-- Modular Arithmetic
--

newtype IntMod (m :: Nat) = IntMod { unwrapN :: Int64 } deriving (Eq)

instance Show (IntMod m) where
  show (IntMod x) = show x

instance KnownNat m => Num (IntMod m) where
  t@(IntMod x) + IntMod y
    | x + y >= modulus = IntMod (x + y - modulus)
    | otherwise = IntMod (x + y)
    where modulus = fromIntegral (natVal t)
  t@(IntMod x) - IntMod y
    | x >= y = IntMod (x - y)
    | otherwise = IntMod (x - y + modulus)
    where modulus = fromIntegral (natVal t)
  t@(IntMod x) * IntMod y = IntMod ((x * y) `rem` modulus)
    where modulus = fromIntegral (natVal t)
  fromInteger n = let result = IntMod (fromInteger (n `mod` fromIntegral modulus))
                      modulus = natVal result
                  in result
  abs = undefined; signum = undefined

{-# RULES
"^9/Int" forall x. x ^ (9 :: Int) = let u = x; v = u * u * u in v * v * v
"^9/Integer" forall x. x ^ (9 :: Integer) = let u = x; v = u * u * u in v * v * v
 #-}

fromIntegral_Int64_IntMod :: KnownNat m => Int64 -> IntMod m
fromIntegral_Int64_IntMod n = result
  where
    result | 0 <= n && n < modulus = IntMod n
           | otherwise = IntMod (n `mod` modulus)
    modulus = fromIntegral (natVal result)

{-# RULES
"fromIntegral/Int->IntMod" fromIntegral = fromIntegral_Int64_IntMod . (fromIntegral :: Int -> Int64) :: Int -> IntMod (10^9 + 7)
"fromIntegral/Int64->IntMod" fromIntegral = fromIntegral_Int64_IntMod :: Int64 -> IntMod (10^9 + 7)
 #-}

instance U.Unboxable (IntMod m) where
  type Rep (IntMod m) = Int64

exEuclid :: (Eq a, Integral a) => a -> a -> (a, a, a)
exEuclid !f !g = loop 1 0 0 1 f g
  where loop !u0 !u1 !v0 !v1 !f 0 = (f, u0, v0)
        loop !u0 !u1 !v0 !v1 !f g =
          case divMod f g of
            (q,r) -> loop u1 (u0 - q * u1) v1 (v0 - q * v1) g r

instance KnownNat m => Fractional (IntMod m) where
  recip t@(IntMod x) = IntMod $ case exEuclid x modulus of
    (1,a,_)  -> a `mod` modulus
    (-1,a,_) -> (-a) `mod` modulus
    _        -> error "not invertible"
    where modulus = fromIntegral (natVal t)
  fromRational = undefined

recipM :: (Eq a, Integral a, Show a) => a -> a -> a
recipM !x modulo = case exEuclid x modulo of
                     (1,a,_) -> a `mod` modulo
                     (-1,a,_) -> (-a) `mod` modulo
                     (g,a,b) -> error $ show x ++ "^(-1) mod " ++ show modulo ++ " failed: gcd=" ++ show g

-- |
-- >>> crt 3 6 2 7
-- 9
-- >>> crt 2 5 3 9
-- 12
crt :: Int64 -> Int64 -> Int64 -> Int64 -> Int64
crt !a1 !m1 !a2 !m2 = let m1' = recipM m1 m2
                          m2' = recipM m2 m1
                      in (m2 * m2' * a1 + m1 * m1' * a2) `mod` (m1 * m2)

max_stars_repo_name: minoki/my-atcoder-solutions

max_issues_repo_name: minoki/my-atcoder-solutions

max_forks_repo_name: minoki/my-atcoder-solutions


text:
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE PatternGuards #-}
{-# LANGUAGE ScopedTypeVariables #-}

module Main
  ( main
  )
where

import Graphics.Gloss hiding (Vector)
import Graphics.Gloss.Interface.Pure.Game hiding (Vector)

import Numeric.LinearAlgebra.Data

import System.Random

import IFS.Base

-- | Main function, when we launch, just start and display the window
main :: IO ()
main = do
  play window white fps initialState render handleEvents stepWorld

-- | Window for the program
window :: Display
window = InWindow "Hacktals" (800, 600) (10, 10)

-- | FPS
fps :: Int
fps = 300

-- | Data structure to hold the state of the fractal
data FractalState = FractalState
  { pointLoc :: Vector R        -- ^ Location of the tracing point
  , points :: [(Float, Float)]  -- ^ Trace of all previous points
  , generator :: StdGen         -- ^ Random generator instance for the "chaos game"
  , autoPlay :: Bool            -- ^ If True, automatically generate the new datapoint
  } deriving Show

-- | Initial state for the program setup
initialState :: FractalState
initialState = FractalState
  { pointLoc = vector [0, 0] -- ^ Start the point at the origin
  , points = [(0, 0)]        -- ^ Starting point
  , generator = mkStdGen 3   -- ^ Start the generator with seed 3
  , autoPlay = False         -- ^ Have autoplay turned off at the start
  }

-- | Rendering function
render :: FractalState -- ^ Fractal info to render
       -> Picture      -- ^ Rendered fractal

render fs = pictures $ map printTrace pts
  where
    pts = points fs
    printTrace (x, y) = translate x y $ color (dark red) $ circleSolid 1

-- | Reactive event handler function
handleEvents :: Event -> FractalState -> FractalState
handleEvents (EventKey (Char 'n') Down _ _) fs = let pl = pointLoc fs
                                                     gen = generator fs
                                                     (pl', (x, y), gen') = chaosStep gen pl
                                                 in fs { pointLoc = pl'
                                                       , generator = gen'
                                                       , points = (x, y):points fs
                                                       }

handleEvents (EventKey (Char 'p') Down _ _) fs = fs { autoPlay = not $ autoPlay fs }
handleEvents _ fs = fs

-- | Stuff to perform at every frame of the animation
stepWorld :: Float -> FractalState -> FractalState
stepWorld _ fs = if autoPlay fs
                 then let pl = pointLoc fs
                          gen = generator fs
                          (pl', (x, y), gen') = chaosStep gen pl
                      in fs { pointLoc = pl'
                            , generator = gen'
                            , points = (x, y):points fs
                            }
                 else fs

chaosStep :: StdGen -> Vector R -> (Vector R, (Float, Float), StdGen)
chaosStep gen pl = let (pl', gen') = chaosIter differentFlake2IFS gen pl
                       (x:y:[]) = toList pl'
                       x' = realToFrac (300 * x)
                       y' = realToFrac (300 * y)
                   in (pl', (x', y'), gen')


text:
{-|
Module: MachineLearning.Clustering
Description: Clustering
Copyright: (c) Alexander Ignatyev, 2017
License: BSD-3
Stability: experimental
Portability: POSIX

Cluster Analysis a.k.a. Clustering - grouping data into coherent subsets.
-}

module MachineLearning.Clustering
(
  Cluster(..)
  , kmeans

  -- * Exported for testing purposes only.
  , kmeansIterM
)

where

import MachineLearning.Types (R, Vector, Matrix)
import Data.List (sortOn, groupBy, minimumBy)
import Control.Applicative ((<$>))
import Control.Monad (forM)
import qualified Control.Monad.Random as RndM
import qualified Data.Vector as V
import qualified Numeric.LinearAlgebra as LA
import MachineLearning.Random (sampleM)

-- | Cluster type (list of samples associtaed with the cluster).
type Cluster = V.Vector Vector

-- | Gets list if the nearest centroid to the sample.
nearestCentroidIndex :: V.Vector Vector  -- ^ list of cluster cetroids;
                     -> Vector            -- ^ sample;
                     -> Int               -- ^ index of the nearest centroid.
nearestCentroidIndex centroids v =
  let distances = V.map (\centroid -> LA.norm_2 (v-centroid)) centroids
  in V.minIndex distances

-- | Calculates cost associated with one claster.
calcClusterCost :: Cluster  -- ^ cluster;
                -> Vector   -- ^ cluster centroid;
                -> R        -- ^ cost value.
calcClusterCost cluster centroid = sum $ fmap (\sample -> LA.norm_2 $ sample-centroid) cluster

-- | Calculates cost value for all clusters.
calcCost :: V.Vector Cluster  -- ^ cluster list;
         -> V.Vector Vector   -- ^ list of cluster centroids;
         -> R                  -- ^ cost value.
calcCost clusters centroids = sum $ V.zipWith calcClusterCost clusters centroids

-- | Calculates centroid (middle point) of the given cluster.
getNewCentroid :: Cluster      -- ^ cluster;
               -> Vector       -- ^ centroid.
getNewCentroid cluster =
  let n = length cluster
      centroid = (sum cluster) / (fromIntegral n)
  in centroid

-- | Calculates new cluster centroids for each cluster.
moveCentroids :: V.Vector Cluster    -- ^ list of clusters;
              -> V.Vector Vector     -- ^ list of cluster centroids.
moveCentroids clusters = fmap getNewCentroid clusters

-- | Build cluster list from list of clusters indices.
buildClusterList :: V.Vector Vector   -- ^ list of samples;
                 -> V.Vector Int      -- ^ list of cluster indices (associated cluster index for each sample);
                 -> V.Vector Cluster  -- ^ list of clusters.
buildClusterList samples clusterIndicesList = V.fromList $ fmap getClusterSamples clusters''
  where clusters' = groupBy (\l r -> snd l == snd r) $ sortOn snd $ zip [0..] $ V.toList clusterIndicesList
        clusters'' = map (map fst) clusters'
        getClusterSamples clusterIndices = V.fromList $ fmap (samples V.!) clusterIndices

-- -- | Run K-Means algorithm once.
kmeansIter :: V.Vector Vector           -- ^ list of samples;
              -> Int                    -- ^ number of clusters (`K`);
              -> V.Vector Vector        -- ^ list of initial centroids;
              -> (V.Vector Cluster, [R])  -- ^ (list of clusters, cost values).
kmeansIter samples k initialCentroids =
  let iter centroids js =
        let clusterIndicesList = fmap (nearestCentroidIndex centroids) samples
            clusters = buildClusterList samples clusterIndicesList
            centroids' = moveCentroids clusters
            j = calcCost clusters centroids'
            diff = sum . fmap LA.norm_2 $ V.zipWith (-) centroids centroids'
        in if diff < 0.001 then (clusters, j:js)
           else iter centroids' (j:js)
  in iter initialCentroids []

-- | Run K-Means algorithm once inside Random Monad.
kmeansIterM :: RndM.RandomGen g =>
               V.Vector Vector  -- ^ list of samples;
               -> Int           -- ^ number of clusters (`K`);
               -> Int           -- ^ iteration number;
               -> RndM.Rand g (V.Vector Cluster, [R])  -- ^ (list of clusters, cost values) inside Random Monad.
kmeansIterM samples k _ = do
  centroids <- sampleM k samples
  return (kmeansIter samples k centroids)

-- | Clusters data using K-Means Algorithm inside Random Monad.
-- Runs K-Means algorithm `N` times, returns the clustering with smaller cost.
kmeans :: RndM.RandomGen g =>
           Int                     -- ^ number of K-Means Algorithm runs (`N`);
           -> Matrix                  -- ^ data to cluster;
           -> Int                     -- ^ desired number of clusters (`K`);
           -> RndM.Rand g (V.Vector Cluster)  -- ^ list of clusters inside Random Monad.
kmeans nIters x k = fst <$>
    (minimumBy (\(_, js1) (_, js2) -> compare (head js1) (head js2))) <$>
    forM [1..nIters] (kmeansIterM samples k)
  where samples = V.fromList $ LA.toRows x

max_stars_repo_name: aligusnet/mltool

max_issues_repo_name: aligusnet/mltool

max_forks_repo_name: aligusnet/mltool


text:
module ComplexLab where

import Data.Complex

-- 1.2 Complex numbers in ~~ Haskell ~~

i :: Num a => a -> Complex a
i n = 0 :+ n

-- Haskell supports complex numbers.
-- The square root of -9, the imaginary number 3i, is written:

squareRootOfNine :: Num a => Complex a
squareRootOfNine =
  i 3 -- Written 3j in python

-- Haskell allows the use of + to add a real number to an imaginary one.
-- We can write the complex solution to `(x − 1)^2 = −9` as:

addRealToImag :: RealFloat a => Complex a
addRealToImag =
  1 + squareRootOfNine -- 1 + 3j in python

-- In fact, the operators +, -, *, /, and ^^ all work with complex numbers.
-- When you add two complex numbers, the real parts are added and the imaginary parts are added.

supportsOperators :: RealFloat a => Complex a
supportsOperators =
  (((1 + i 3) / 8) ^^ 2) * i 3
  -- ((((1 + 3j) / 8) ** 2) * 3j) in python

max_stars_repo_name: danbroooks/coding-the-matrix

max_issues_repo_name: danbroooks/coding-the-matrix

max_forks_repo_name: danbroooks/coding-the-matrix


text:
{-# LANGUAGE ViewPatterns #-}
--------------------------------------------------------------------------------
-- |
-- Module      : ArrayFire.Signal
-- Copyright   : David Johnson (c) 2019-2020
-- License     : BSD 3
-- Maintainer  : David Johnson <djohnson.m@gmail.com>
-- Stability   : Experimental
-- Portability : GHC
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft.htm)
--
-- Functions performing signal processing on 'Array'
--
-- >>> input = vector 3 [10,20,30]
-- >>> positions = vector 5 [0.0, 0.5, 1.0, 1.5, 2.0]
-- >>> approx1 @Double input positions Cubic 0.0
-- ArrayFire Array
-- [5 1 1 1]
--    10.0000
--    13.7500
--    20.0000
--    26.2500
--    30.0000
--
--------------------------------------------------------------------------------
module ArrayFire.Signal where

import Data.Complex

import ArrayFire.FFI
import ArrayFire.Internal.Signal
import ArrayFire.Internal.Types

-- | 'approx1' interpolates data along the first dimensions
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__approx1.htm)
--
-- Interpolation is performed assuming input data is equally spaced with indices in the range [0, n). The positions are sampled with respect to data at these locations.
--
-- >>> input = vector 3 [10,20,30]
-- >>> positions = vector 5 [0.0, 0.5, 1.0, 1.5, 2.0]
-- >>> approx1 @Double input positions Cubic 0.0
-- ArrayFire Array
-- [5 1 1 1]
--    10.0000
--    13.7500
--    20.0000
--    26.2500
--    30.0000
approx1
  :: AFType a
  => Array a
  -- ^ the input array
  -> Array a
  -- ^ array contains the interpolation locations
  -> InterpType
  -- ^ is the interpolation type, it can take one of the values defined by 'InterpType'
  -> Float
  -- ^ is the value that will set in the output array when certain index is out of bounds
  -> Array a
  -- ^ is the array with interpolated values
approx1 arr1 arr2 (fromInterpType -> i1) f =
  op2 arr1 arr2 (\p x y -> af_approx1 p x y i1 f)

-- | approx2 performs interpolation on data along the first and second dimensions.
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__approx2.htm)
--
-- Interpolation is performed assuming input data is equally spaced with indices in the range [0, n) along each dimension. The positions are sampled with respect to data at these locations.
--
-- >>> input = matrix @Double (3,3) [ [ 1.0,1.0,1.0 ], [ 2.0, 2.0, 2.0 ], [ 3.0,3.0,3.0 ] ]
-- >>> positions1 = matrix @Double (2,2) [ [ 0.5,1.5 ],[ 0.5,1.5 ] ]
-- >>> positions2 = matrix @Double (2,2) [ [ 0.5,0.5 ],[ 1.5,1.5 ] ]
-- >>> approx2 @Double input positions1 positions2 Cubic 0.0
-- ArrayFire Array
-- [2 2 1 1]
--     1.3750     2.6250
--     1.3750     2.6250
--
approx2
  :: AFType a
  => Array a
  -- ^ is the input array
  -> Array a
  -- ^ array contains the interpolation locations for first dimension
  -> Array a
  -- ^ array contains the interpolation locations for second dimension
  -> InterpType
  -- ^ is the interpolation type, it can take one of the values defined by 'InterpType'
  -> Float
  -- ^ is the value that will set in the output array when certain index is out of bounds
  -> Array a
  -- ^	is the array with interpolated values
approx2 arr1 arr2 arr3 (fromInterpType -> i1) f =
  op3 arr1 arr2 arr3 (\p x y z -> af_approx2 p x y z i1 f)

-- DMJ: Where did these functions go? Were they removed?
-- http://arrayfire.org/docs/group__approx__mat.htm
-- approx1Uniform
--   :: AFType a
--   => Array a
--   -> Array a
--   -> Int
--   -> Double
--   -> Double
--   -> InterpType
--   -> Float
--   -> Array a
-- approx1Uniform arr1 arr2 (fromIntegral -> i1) d1 d2 (fromInterpType -> interp) f =
--   op2 arr1 arr2 (\p x y -> af_approx1_uniform p x y i1 d1 d2 interp f)

-- approx2Uniform
--   :: AFType a
--   => Array a
--   -> Array a
--   -> Int
--   -> Double
--   -> Double
--   -> Array a
--   -> Int
--   -> Double
--   -> Double
--   -> InterpType
--   -> Float
--   -> Array a
-- approx2Uniform arr1 arr2 (fromIntegral -> i1) d1 d2 arr3 (fromIntegral -> i2) d3 d4 (fromInterpType -> interp) f =
--   op3 arr1 arr2 arr3 (\p x y z -> af_approx2_uniform p x y i1 d1 d2 z i2 d3 d4 interp f)

-- | Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft.htm#ga64d0db9e59c9410ba738591ad146a884)
--
-- The Fast Fourier Transform (FFT) is an efficient algorithm to compute the discrete Fourier transform (DFT) of a signal or array. This is most commonly used to convert data in the time (or space) domain to the frequency domain, Then, the inverse FFT (iFFT) is used to return the data to the original domain.
--
-- >>> fft (vector @Double 10 [1..]) 2.0 10
-- ArrayFire Array
-- [10 1 1 1]
--          (110.0000,0.0000)
--          (-10.0000,30.7768)
--          (-10.0000,13.7638)
--          (-10.0000,7.2654)
--          (-10.0000,3.2492)
--          (-10.0000,0.0000)
--          (-10.0000,-3.2492)
--          (-10.0000,-7.2654)
--          (-10.0000,-13.7638)
--          (-10.0000,-30.7768)
fft
  :: (AFType a, Fractional a)
  => Array a
  -- ^ input 'Array'
  -> Double
  -- ^ the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals - used to either truncate or pad the input signals.
  -> Array a
  -- ^ is the transformed array
fft a d (fromIntegral -> x) =
  op1 a (\j k -> af_fft j k d x)

-- | Fast Fourier Transform (in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft.htm#gaa2f03c9ee1cb80dc184c0b0a13176da1)
--
-- C Interface for fast fourier transform on one dimensional signals.
--
-- *Note* The input in must be a complex array
--
fftInPlace
  :: (AFType a, Fractional a)
  => Array (Complex a)
  -- ^ is the input array on entry and the output of 1D forward fourier transform at exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
fftInPlace a d = a `inPlace` (flip af_fft_inplace d)

-- | Fast Fourier Transform (2-dimensional)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft2.htm#gaab3fb1ed398e208a615036b4496da611)
--
-- C Interface for fast fourier transform on two dimensional signals.
--
fft2
  :: AFType a
  => Array a
  -- ^ the input array
  -> Double
  -- ^ the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Array a
  -- ^ the transformed array
fft2 a d x y =
  op1 a (\j k -> af_fft2 j k d (fromIntegral x) (fromIntegral y))

-- | Fast Fourier Transform (2-dimensional, in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft2.htm#gacdeebb3f221ae698833dc4900a172b8c)
--
-- C Interface for fast fourier transform on two dimensional signals.
--
-- *Note* The input in must be a complex array
--
fft2_inplace
  :: (Fractional a, AFType a)
  => Array (Complex a)
  -- ^  input array on entry and the output of 2D forward fourier transform on exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
fft2_inplace a d = a `inPlace` (flip af_fft2_inplace d)

-- | Fast Fourier Transform (3-dimensional)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft3.htm#ga5138ef1740ece0fde2c796904d733c12)
--
-- C Interface for fast fourier transform on three dimensional signals.
--
fft3
  :: AFType a
  => Array a
  -- ^ the input array
  -> Double
  -- ^ the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along third dimension - used to either truncate/pad the input
  -> Array a
  -- ^ the transformed array
fft3 a d x y z =
  op1 a (\j k -> af_fft3 j k d (fromIntegral x) (fromIntegral y) (fromIntegral z))

-- | Fast Fourier Transform (3-dimensional, in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft2.htm#gacdeebb3f221ae698833dc4900a172b8c)
--
-- C Interface for fast fourier transform on three dimensional signals.
--
-- *Note* The input in must be a complex array
--
fft3_inplace
  :: (Fractional a, AFType a)
  => Array (Complex a)
  -- ^  input array on entry and the output of 3D forward fourier transform on exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
fft3_inplace a d = a `inPlace` (flip af_fft3_inplace d)

-- | Inverse Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft.htm#ga2d62c120b474b3b937b0425c994645fe)
--
-- C Interface for inverse fast fourier transform on one dimensional signals.
--
ifft
  :: AFType a
  => Array a
  -- ^ the input array
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^  is the length of output signals - used to either truncate or pad the input signals
  -> Array a
  -- ^ the transformed array
ifft a d x =
  op1 a (\j k -> af_ifft j k d (fromIntegral x))

-- | Inverse Fast Fourier Transform (in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft.htm#ga827379bef0e2cadb382c1b6301c91429)
--
-- C Interface for fast fourier transform on one dimensional signals.
--
ifft_inplace
  :: (AFType a, Fractional a)
  => Array (Complex a)
  -- ^ is the input array on entry and the output of 1D forward fourier transform at exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
ifft_inplace a d = a `inPlace` (flip af_ifft_inplace d)

-- | Inverse Fast Fourier Transform (2-dimensional signals)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft2.htm#ga7cd29c6a35c19240635b62cc5c30dc4f)
--
-- C Interface for inverse fast fourier transform on two dimensional signals.
--
ifft2
  :: AFType a
  => Array a
  -- ^ the input array
  -> Double
  -- ^ the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Array a
  -- ^ the transformed array
ifft2 a d x y =
  op1 a (\j k -> af_ifft2 j k d (fromIntegral x) (fromIntegral y))

-- | Inverse Fast Fourier Transform (2-dimensional, in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft2.htm#ga9e6a165d44306db4552a56d421ce5d05)
--
-- C Interface for fast fourier transform on two dimensional signals.
--
ifft2_inplace
  :: (AFType a, Fractional a)
  => Array (Complex a)
  -- ^ is the input array on entry and the output of 1D forward fourier transform at exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
ifft2_inplace a d = a `inPlace` (flip af_ifft2_inplace d)

-- | Inverse Fast Fourier Transform (3-dimensional)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft3.htm)
--
-- C Interface for inverse fast fourier transform on three dimensional signals.
--
ifft3
  :: AFType a
  => Array a
  -- ^ the input array
  -> Double
  -- ^ the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along third dimension - used to either truncate/pad the input
  -> Array a
  -- ^ the transformed array
ifft3 a d x y z =
  op1 a (\j k -> af_ifft3 j k d (fromIntegral x) (fromIntegral y) (fromIntegral z))

-- | Inverse Fast Fourier Transform (3-dimensional, in-place)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__ifft3.htm#ga439a7a49723bc6cf77cf4fe7f8dfe334)
--
-- C Interface for fast fourier transform on two dimensional signals.
--
ifft3_inplace
  :: (AFType a, Fractional a)
  => Array (Complex a)
  -- ^ is the input array on entry and the output of 1D forward fourier transform at exit
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> IO ()
ifft3_inplace a d = a `inPlace` (flip af_ifft3_inplace d)

-- | Real to Complex Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__r2c.htm#ga7486f342182a18e773f14cc2ab4cb551)
--
-- C Interface for real to complex fast fourier transform for one dimensional signals.
--
-- The first dimension of the output will be of size (pad0 / 2) + 1. The second dimension of the output will be pad1. The third dimension of the output will be pad 2.
--
fftr2c
  :: (Fractional a, AFType a)
  => Array a
  -- ^ is a real array
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Array a
  -- ^ is a complex array containing the non redundant parts of in.
fftr2c a d x =
  op1 a (\j k -> af_fft_r2c j k d (fromIntegral x))

-- | Real to Complex Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__r2c.htm#ga7486f342182a18e773f14cc2ab4cb551)
--
-- C Interface for real to complex fast fourier transform for two dimensional signals.
--
-- The first dimension of the output will be of size (pad0 / 2) + 1. The second dimension of the output will be pad1. The third dimension of the output will be pad 2.
--
fft2r2c
  :: (Fractional a, AFType a)
  => Array a
  -- ^ is a real array
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Array a
  -- ^ is a complex array containing the non redundant parts of in.
fft2r2c a d x y =
  op1 a (\j k -> af_fft2_r2c j k d (fromIntegral x) (fromIntegral y))

-- | Real to Complex Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__r2c.htm#gab4ca074b54218b74d8cfbda63d38be51)
--
-- C Interface for real to complex fast fourier transform for three dimensional signals.
--
-- The first dimension of the output will be of size (pad0 / 2) + 1. The second dimension of the output will be pad1. The third dimension of the output will be pad 2.
--
fft3r2c
  :: (Fractional a, AFType a)
  => Array a
  -- ^ is a real array
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Int
  -- ^ is the length of output signals along first dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along second dimension - used to either truncate/pad the input
  -> Int
  -- ^ is the length of output signals along third dimension - used to either truncate/pad the input
  -> Array a
  -- ^ is a complex array containing the non redundant parts of in.
fft3r2c a d x y z =
  op1 a (\j k -> af_fft3_r2c j k d (fromIntegral x) (fromIntegral y) (fromIntegral z))

-- | Complex to Real Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__c2r.htm#gaa5efdfd84213a4a07d81a5d534cde5ac)
--
-- C Interface for complex to real fast fourier transform for one dimensional signals.
--
-- The first dimension of the output will be 2 * dim0 - 1 if is_odd is true else 2 * dim0 - 2 where dim0 is the first dimension of the input. The remaining dimensions are unchanged.
--
fftc2r
  :: AFType a
  => Array a
  -- ^ is a complex array containing only the non redundant parts of the signals.
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Bool
  -- ^ is a flag signifying if the output should be even or odd size
  -> Array a
  -- ^ is a real array containing the output of the transform.
fftc2r a cm (fromIntegral . fromEnum -> cd) = op1 a (\x y -> af_fft_c2r x y cm cd)

-- | Complex to Real Fast Fourier Transform (2-dimensional)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__c2r.htm#gaaa7da16f226cacaffced631e08da4493)
--
-- C Interface for complex to real fast fourier transform for two dimensional signals.
--
-- The first dimension of the output will be 2 * dim0 - 1 if is_odd is true else 2 * dim0 - 2 where dim0 is the first dimension of the input. The remaining dimensions are unchanged.
--
fft2C2r
  :: AFType a
  => Array a
  -- ^ is a complex array containing only the non redundant parts of the signals.
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Bool
  -- ^ is a flag signifying if the output should be even or odd size
  -> Array a
  -- ^ is a real array containing the output of the transform.
fft2C2r a cm (fromIntegral . fromEnum -> cd) = op1 a (\x y -> af_fft2_c2r x y cm cd)

-- | Complex to Real Fast Fourier Transform (3-dimensional)
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__c2r.htm#gaa9b3322d9ffab15268919e1f114bed24)
--
-- C Interface for complex to real fast fourier transform for three dimensional signals.
--
-- The first dimension of the output will be 2 * dim0 - 1 if is_odd is true else 2 * dim0 - 2 where dim0 is the first dimension of the input. The remaining dimensions are unchanged.
--
fft3C2r
  :: AFType a
  => Array a
  -- ^ is a complex array containing only the non redundant parts of the signals.
  -> Double
  -- ^ is the normalization factor with which the input is scaled after the transformation is applied
  -> Bool
  -- ^ is a flag signifying if the output should be even or odd size
  -> Array a
  -- ^ is a real array containing the output of the transform.
fft3C2r a cm (fromIntegral . fromEnum -> cd) = op1 a (\x y -> af_fft3_c2r x y cm cd)

-- | Convolution Integral for one dimensional data
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__convolve1.htm#ga25d77b794463b5cd72cd0b7f4af140d7)
--
-- C Interface for convolution on one dimensional signals.
--
-- *Note* The default parameter of domain, AF_CONV_AUTO, heuristically switches between frequency and spatial domain.
--
convolve1
  :: AFType a
  => Array a
  -- ^ the input signal
  -> Array a
  -- ^ the signal that shall be flipped for the convolution operation
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> ConvDomain
  -- ^ specifies if the convolution should be performed in frequency os spatial domain
  -> Array a
  -- ^ convolved array
convolve1 a b (toConvMode -> cm) (fromConvDomain -> cd) = op2 a b (\x y z -> af_convolve1 x y z cm cd)

-- | Convolution Integral for two dimensional data
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__convolve2.htm#ga25d77b794463b5cd72cd0b7f4af140d7)
--
-- C Interface for convolution on two dimensional signals.
--
-- *Note* The default parameter of domain, AF_CONV_AUTO, heuristically switches between frequency and spatial domain.
--
convolve2
  :: AFType a
  => Array a
  -- ^ the input signal
  -> Array a
  -- ^ the signal that shall be flipped for the convolution operation
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> ConvDomain
  -- ^ specifies if the convolution should be performed in frequency os spatial domain
  -> Array a
  -- ^ convolved array
convolve2 a b (toConvMode -> cm) (fromConvDomain -> cd) = op2 a b (\x y z -> af_convolve2 x y z cm cd)

-- | Convolution Integral for three dimensional data
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__convolve3.htm#ga25d77b794463b5cd72cd0b7f4af140d7)
--
-- C Interface for convolution on three dimensional signals.
--
-- *Note* The default parameter of domain, AF_CONV_AUTO, heuristically switches between frequency and spatial domain.
--
convolve3
  :: AFType a
  => Array a
  -- ^ the input signal
  -> Array a
  -- ^ the signal that shall be flipped for the convolution operation
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> ConvDomain
  -- ^ specifies if the convolution should be performed in frequency os spatial domain
  -> Array a
  -- ^ convolved array
convolve3 a b (toConvMode -> cm) (fromConvDomain -> cd) =
  op2 a b (\x y z -> af_convolve3 x y z cm cd)

-- | C Interface for separable convolution on two dimensional signals.
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__convolve.htm#gaeb6ba88155cf3ef29d93f97b147e372f)
--
-- C Interface for separable convolution on two dimensional signals.
--
-- *Note* The default parameter of domain, AF_CONV_AUTO, heuristically switches between frequency and spatial domain.
--
convolve2Sep
  :: AFType a
  => Array a
  -- ^ filter that has to be applied along the coloumns
  -> Array a
  -- ^ filter that has to be applied along the rows
  -> Array a
  -- ^ the input array
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> Array a
  -- ^ convolved array
convolve2Sep a b c (toConvMode -> d) = op3 a b c (\x y z j -> af_convolve2_sep x y z j d)

-- DMJ: did this get removed? Can't find in latest docs
-- fftConvolve1
--   :: AFType a
--   => Array a
--   -> Array a
--   -> ConvMode
--   -> Array a
-- fftConvolve1 a b (toConvMode -> c) = op2 a b (\x y z -> af_fft_convolve1 x y z c)

-- | 2D Convolution using Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__convolve2.htm#gab52ebe631d8358cdef1b5c8a95550556)
--
-- A convolution is a common operation between a source array, a, and a filter (or kernel) array b. The answer to the convolution is the same as computing the coefficients in polynomial multiplication, if a and b are the coefficients.
--
-- C Interface for FFT-based convolution on two dimensional signals.
--
fftConvolve2
  :: AFType a
  => Array a
  -- ^ is the input signal
  -> Array a
  --  ^ is the signal that shall be used for the convolution operation
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> Array a
  -- ^  is convolved array
fftConvolve2 a b (toConvMode -> c) = op2 a b (\x y z -> af_fft_convolve2 x y z c)

-- | 3D Convolution using Fast Fourier Transform
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft__convolve3.htm)
--
-- A convolution is a common operation between a source array, a, and a filter (or kernel) array b. The answer to the convolution is the same as computing the coefficients in polynomial multiplication, if a and b are the coefficients.
--
-- C Interface for FFT-based convolution on three dimensional signals.
--
fftConvolve3
  :: AFType a
  => Array a
  -- ^ is the input signal
  -> Array a
  --  ^ is the signal that shall be used for the convolution operation
  -> ConvMode
  -- ^ indicates if the convolution should be expanded or not(where output size equals input)
  -> Array a
  -- ^  is convolved array
fftConvolve3 a b (toConvMode -> c) = op2 a b (\x y z -> af_fft_convolve3 x y z c)

-- | Finite Impulse Filter.
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fir.htm#ga2a850e69775eede4709e0d607bba240b)
--
-- Finite impulse filters take an input x and a co-efficient array b to generate an output y such that:
--
-- C Interface for finite impulse response filter.
--
fir
  :: AFType a
  => Array a
  -- ^ is the input signal to the filter
  -> Array a
  -- ^ is the array containing the coefficients of the filter
  -> Array a
  -- ^ is the output signal from the filter
fir a b = op2 a b af_fir

-- | Infinite Impulse Filter.
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__iir.htm#ga7adcc364da0a66cdfd2bb351215456c4)
--
-- Infinite impulse filters take an input x and a feedforward array b, feedback array a to generate an output y such that:
--
-- C Interface for infinite impulse response filter.
--
-- *Note* The feedforward coefficients are currently limited to a length of 512
--
iir
  :: AFType a
  => Array a
  -- ^ the array containing the feedforward coefficient
  -> Array a
  -- ^ is the array containing the feedback coefficients
  -> Array a
  -- ^ is the input signal to the filter
  -> Array a
  -- ^ the output signal from the filter
iir a b c = op3 a b c af_iir

-- | Median Filter
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__image__func__medfilt.htm#gaaf3f62f2de0f4dc315b831e494e1b2c0)
--
-- A median filter is similar to the arbitrary filter except that instead of a weighted sum, the median value of the pixels covered by the kernel is returned.
--
-- C Interface for median filter.
--
medFilt
  :: AFType a
  => Array a
  -- ^ 'Array' is the input image
  -> Int
  -> Int
  -> BorderType
  -> Array a
  -- ^ 'Array' is the processed image
medFilt a l w (fromBorderType -> b) =
 a `op1` (\x y -> af_medfilt x y (fromIntegral l) (fromIntegral w) b)

-- | 1D Median Filter
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__image__func__medfilt.htm#gad108ea62cbbb5371bd14a17d06384359)
--
-- A median filter is similar to the arbitrary filter except that instead of a weighted sum, the median value of the pixels covered by the kernel is returned.
--
-- C Interface for 1D median filter.
--
medFilt1
  :: AFType a
  => Array a
  -- ^ 'Array' is the input signal
  -> Int
  -- ^ Is the kernel width
  -> BorderType
  -- ^ value will decide what happens to border when running filter in their neighborhood. It takes one of the values [AF_PAD_ZERO | AF_PAD_SYM]
  -> Array a
  -- ^ 'Array' is the processed signal
medFilt1 a w (fromBorderType -> b) =
 a `op1` (\x y -> af_medfilt1 x y (fromIntegral w) b)

-- | 2D Median Filter
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__image__func__medfilt.htm#ga2cb99dca5842f74f6b9cd28eb187a9cd)
--
-- A median filter is similar to the arbitrary filter except that instead of a weighted sum, the median value of the pixels covered by the kernel is returned.
--
-- C Interface for 2D median filter.
--
medFilt2
  :: AFType a
  => Array a
  -- ^ 'Array' is the input image
  -> Int
  -- ^ the kernel height
  -> Int
  -- ^ the kernel width
  -> BorderType
  -- ^ value will decide what happens to border when running filter in their neighborhood. It takes one of the values [AF_PAD_ZERO | AF_PAD_SYM]
  -> Array a
  -- ^ 'Array' is the processed image
medFilt2 a l w (fromBorderType -> b) =
 a `op1` (\x y -> af_medfilt2 x y (fromIntegral l) (fromIntegral w) b)

-- | C Interface for setting plan cache size.
--
-- [ArrayFire Docs](http://arrayfire.org/docs/group__signal__func__fft.htm#ga4ddef19b43d9a50c97b1a835df60279a)
--
-- This function doesn't do anything if called when CPU backend is active. The plans associated with the most recently used array sizes are cached.
--
-- >>> setFFTPlanCacheSize 2
-- ()
--
setFFTPlanCacheSize
  :: Int
  -- ^ is the number of plans that shall be cached
  -> IO ()
setFFTPlanCacheSize =
 afCall . af_set_fft_plan_cache_size . fromIntegral

max_stars_repo_name: arrayfire/arrayfire-haskell

max_issues_repo_name: arrayfire/arrayfire-haskell

max_forks_repo_name: arrayfire/arrayfire-haskell