text:
{-# LANGUAGE DuplicateRecordFields     #-}
{-# LANGUAGE ExistentialQuantification #-}
{-# LANGUAGE RankNTypes                #-}
{-# LANGUAGE ScopedTypeVariables       #-}
{-# LANGUAGE TemplateHaskell           #-}

module Minimizer(module Minimizer) where

import           Brain
import           Control.Arrow
import           Control.Lens
import           Control.Monad.Random
import           Convenience
import qualified Data.Vector.Storable           as V
import           GHC.TypeLits
import           Numeric.GSL.Minimization
import           Numeric.GSL.SimulatedAnnealing
import           Simulator
import qualified Graphics.Rendering.Chart.Easy  as E
import qualified Graphics.Rendering.Chart.Gtk   as C
import           Numeric.LinearAlgebra.Data
import           System.IO.Unsafe
import qualified Numeric.AD as AD

type Minimizer n = (forall a. Floating a => Weights n a -> a) -> Weights n Double -> Weights n Double

newtype NelderMeadsParams = NelderMeadsParams Int

instance Default NelderMeadsParams where
  auto = NelderMeadsParams 1000

instance Default SimulatedAnnealingParams where
  auto = SimulatedAnnealingParams 100 100 1.0 1.0 100 1.5 10

newtype AutoDiffParams = AutoDiffParams Int

instance Default AutoDiffParams where
  auto = AutoDiffParams 10

data MinSettings = Annealing SimulatedAnnealingParams | NelderMeads NelderMeadsParams | AutoDiff AutoDiffParams

instance Default MinSettings where
  auto = AutoDiff auto

data Settings = Settings {
    _minSettings      :: MinSettings,
    _seed             :: Int,
    _systems          :: Int,
    _groupedBy        :: Int,
    _iterRange        :: (Int,Int),
    _numMinimizations :: Int
}
makeLenses ''Settings

instance Default Settings where
  auto = Settings {
    _minSettings = auto,
    _seed = 4374653543,
    _systems = 3,
    _groupedBy = 30,
    _iterRange = (400,400),
    _numMinimizations = 5
  }

data NeuralSim system enabled all f = NeuralSim {
    _settings   :: Settings,
    _weights    :: Weights enabled f,
    _restorer   :: Restorer enabled all f,
    _setWeights :: forall a . Floating a => Weights enabled a -> system a -> system a
}
makeLenses ''NeuralSim

minimizeS :: (KnownNat n) => NelderMeadsParams -> Minimizer n
minimizeS (NelderMeadsParams iterations) cost xi =
  minimizeV NMSimplex2 0.0001 iterations (V.replicate (ssize xi) 0.2) (fromVec &. cost) (toVec xi)
  -- & chartify
  & fst & fromVec

chartify :: (a, Matrix Double) -> (a, ())
chartify = second (toColumns &. (!!1)
            &. Numeric.LinearAlgebra.Data.toList
            &. zip [1::Int ..]
            &. (:[]) &. E.line "30" &. E.plot
            &. C.toWindow 500 500 &. unsafePerformIO)
            &. (\t -> seq (snd t) t)

annealing :: (KnownNat n) => SimulatedAnnealingParams -> Minimizer n
annealing anParams cost xi = simanSolve 123 (ssize xi) anParams xi cost metricDist stepFunction (Just $ const "")
  where metricDist a1 a2 = sZipWith (\x1 x2 -> (x1-x2)*(x1-x2)) a1 a2 & ssum & sqrt
        stepFunction rands stepSize current = rands & fromVec &> (\x -> x*2*stepSize - stepSize) & sZipWith (+) current

automaticDiff :: (KnownNat n) => AutoDiffParams -> Minimizer n
automaticDiff (AutoDiffParams iterations) cost (Sized weights) = AD.conjugateGradientDescent (Sized &. cost) weights & take iterations & last & Sized

toMinimizer :: KnownNat n => MinSettings -> Minimizer n
toMinimizer (NelderMeads pars) = minimizeS pars
toMinimizer (Annealing pars)   = annealing pars
toMinimizer (AutoDiff pars)    = automaticDiff pars

minimizer :: (KnownNat all,KnownNat enabled,Simulator system) => NeuralSim system enabled all Double -> IO ()
minimizer n = print (evalRand (minimizeRand n) (mkStdGen $ n^.settings.seed))

minimizeRand :: (KnownNat all, KnownNat enabled, Simulator system) =>
                    NeuralSim system enabled all Double -> Rand StdGen ([(Double,Double)], Weights all Double)
minimizeRand n = composeN (n^.settings.numMinimizations) episodic ([], n^.weights)
                 &> second (n^.restorer)
  where
    episodic (costs,weights) =
      episodicM (toMinimizer (n^.settings.minSettings)) (n^.settings) (n^.setWeights) weights
      &> first (:costs)

episodicM :: forall s enabled. (Simulator s)
          => Minimizer enabled
          -> Settings
          -> (forall a . (Floating a) => Weights enabled a -> s a -> s a)
          -> Weights enabled Double
          -> Rand StdGen ((Double,Double), Weights enabled Double)
episodicM optimizer settins weightSetter initialWeights = do
  iters <- getRandomR (settins^.iterRange)
  randSystems <- getRandoms &> take (settins^.systems)

  let episodeSize = settins^.groupedBy

  let simsEpisode :: Floating a => [s a] -> [s a]
      simsEpisode = apply episodeSize simsStep

  let episodes = iters `div` episodeSize

  let episodedSystems :: Floating a => [[s a]]
      episodedSystems = randSystems
                        &> realToFracSim
                        &  iterate simsEpisode
                        &  take episodes

  let costOfWeights :: (Floating a) => Weights enabled a -> a
      costOfWeights ws = episodedSystems
                         &>>weightSetter ws
                         &> simsEpisode
                         &> simsCost
                         &  sum

  let resultWeights = optimizer costOfWeights initialWeights
  return ((costOfWeights initialWeights,costOfWeights resultWeights), resultWeights)

-- inputRecreator :: _ => Brain n n numWeights -> Minimizer numWeights -> g (Weights numWeights)
-- inputRecreator (Brain feed) optimizer =
--   getRandoms &>> randWeights
--              &> (\inputs -> optimizer (\w -> let outputs = inputs &> feed w
--                                              in zipWith hamming inputs outputs & sum)
--                                       (randWeights 23423))
--   where
--     hamming :: Weights n -> Weights n -> Double
--     hamming w1 w2 = sZipWith (\x y -> (x-y)**2) w1 w2 & sum

max_stars_repo_name: bmabsout/neural-swarm

max_issues_repo_name: bmabsout/neural-swarm

max_forks_repo_name: bmabsout/neural-swarm


text:
module Main where

import Options.Applicative

import Pipes
import qualified Pipes.Prelude as P

import System.Random.MWC

import Statistics.Quantile.Bench.Gen (sampleGaussian)

data Params = Params Int Double Double
  deriving (Eq, Show)

params = Params <$> argument auto (metavar "COUNT")
                <*> argument auto (metavar "MEAN")
                <*> argument auto (metavar "STDDEV")

args = info (helper <*> params)
            (   fullDesc
             <> progDesc "Write COUNT real numbers drawn from a Gaussian."
            )

main :: IO ()
main = do
    ps <- execParser args
    gen <- createSystemRandom
    runEffect $ gauss gen ps >-> P.map show >-> P.stdoutLn
  where
    gauss gen (Params n mu sigma) = P.replicateM n $ sampleGaussian gen mu sigma

ext: hs

max_stars_repo_name: fractalcat/slides

max_issues_repo_name: fractalcat/slides

max_forks_repo_name: fractalcat/slides


text:
{-
#!/usr/bin/env stack
--stack --resolver lts-11.3 --install-ghc runghc --stack-yaml /home/aku/kanazawa/dev/cpvoh/stack.yaml
-}

{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE Strict #-}
{-# LANGUAGE FlexibleContexts #-}

--import CPVO.IO.Plot.Band

import CPVO.IO -- inshell2text
--import System.Environment (getArgs)
import Data.List.Split (splitOn)
import qualified Data.Text as T
import qualified Data.Text.Read as T
import Data.List -- intercalate
import Data.Either -- rights
import Numeric.LinearAlgebra.Data -- toLists, fromLists
import Data.Maybe -- fromJust
import Data.String  -- IsString
import System.Directory
import Options.Applicative as O
import Data.Semigroup -- <>

main :: IO ()
main = do
  (Opts fOut useOldBw judulUtama yr atomOs daftarLengkap) <- execParser withHelp
  plotBand fOut useOldBw judulUtama yr atomOs
    $ map (splitOn "#") daftarLengkap
  putStrLn "========DONE======="

debugLine :: Show a => Bool -> a -> IO ()
debugLine False _ = return ()
debugLine _ ss = putStrLn $ "===debug=== " ++ show ss

plotInit :: String
plotInit = unlines $ "set datafile missing '-'":
                     "plot 0 lt -1 lc rgb 'black' title '' , \\":[]

endMultiplot :: String
endMultiplot = unlines [ "unset multiplot" ]

plotBand :: String -> Bool -> String -> String -> String -> [[String]]
         -> IO ()
plotBand fOut useOldBw judulUtama yr atomOs daftarLengkap = do
  tempDir <- (T.unpack . head) <$> inshell2text "mktemp -d -p ./"

  putStrLn $ unlines $ concat daftarLengkap
  (xrangeatas,ticksbaru,arrow,plotplate) <- genSinglePic useOldBw atomOs daftarLengkap ("","",Nothing,[])

  putStrLn $ "============TEMPGLT"
  tempGLT <- head <$> inshell2text "mktemp -p ./"
  putStrLn $ show tempGLT

  writeFile (T.unpack tempGLT) $ genTEMPGLT tempDir judulUtama yr xrangeatas ticksbaru (fromMaybe "" arrow) plotplate

  system_ $ "gnuplot " ++ T.unpack tempGLT
  let target = map (\x -> unwords [ "mv "
                                  , "hasil" ++ x
                                  , "../" ++ fOut ++ x
                                  ]) [".eps",".png"]
  withCurrentDirectory tempDir $
    mapM_ system_ $
      "ps2eps --rotate=+ tmp.ps":
      "epstool --copy -b --quiet tmp.eps hasil.eps":
      "epstopdf hasil.eps":
      "pdftocairo -r 150 -singlefile -jpeg hasil.pdf tmp":
      "convert tmp.jpg -rotate 0 hasil.png":
      "rm -f tmp.jpg":
      target
  putStrLn "===end  : plotBand===="

genSinglePic :: Bool -> String -> [[String]]
              ->    (String, String, Maybe String, [String])
              -> IO (String, String, Maybe String, [String])
genSinglePic _        _      []            res              = return res
genSinglePic useOldBw atomOs (daftarLengkap:ds) (_,_,_,res) = do
  allPBAND@((xrangeatas,ticksbaru,arrow,_,_):_) <- plotSinglePic daftarLengkap useOldBw atomOs 1 []
  let generatedPFBAND = (map takePFBAND allPBAND)
  let generatedPBAND = T.intercalate "," $ filter (not . T.null) $ map fst generatedPFBAND
  let generatedFATBAND = T.intercalate "," $ filter (not . T.null) $ map snd generatedPFBAND
  genSinglePic useOldBw atomOs ds (xrangeatas, ticksbaru, arrow, (T.unpack $ T.intercalate "," $ filter (not . T.null) [generatedPBAND, generatedFATBAND]):res)
    where
      takePFBAND (_,_,_,p,f) = (p,f)

plotSinglePic :: [String]
               -> Bool
               -> String
               -> Int
               -> [(String, String, Maybe String, T.Text, T.Text)]
               -> IO [(String, String, Maybe String, T.Text, T.Text)]

plotSinglePic []                   _        _       _      res = return res
plotSinglePic (daftarLengkap:sisa) useOldBw atomOs colorId res = do
  let (foldernya:spinnya:legend:_) = splitOn ":" daftarLengkap
  bandFiles <- inshell2text $ unwords ["ls",concat[foldernya,"/bnd*spin",spinnya]]
  let daftarFolder = ":" ++ intercalate "@" [legend,spinnya,foldernya]

  let generatedPBAND = T.intercalate ", " $ map (genPlotPlate colorId) $ map (\x -> T.concat["'",x,"'"]) bandFiles
  valBand@(valBandX:valBandY:_) <- fmap (words . T.unpack . head) $ inshell2text $ unwords ["cat",concat [foldernya,"/bnd*spin*",spinnya],"|sed -e '/^#/d' |awk '{if ($3>0.1) print $2,$3}'|sort -k 2n|head -1"]
  (_:condBandY:_) <- fmap (words . T.unpack . head) $ inshell2text $ unwords ["cat", concat[foldernya,"/bnd*spin*",spinnya],"| sed -e '/^#/d' |awk '{if ($3<=0) print $2,$3}'|sort -k 2nr -u|sed -e '/^ *$/d'|head -1"]
  putStrLn $ "===valBand==" ++ show (valBand :: [String])
  gapCoordBandGap <- fmap runGAPband $ inshell2text $ "cat " ++ concat [foldernya,"/bnd*spin*",spinnya]
  let bandGap = last $ words gapCoordBandGap
  putStrLn $ "===allBand=== " ++ bandGap
  let arrow = genArrow bandGap valBandX valBandY condBandY
--  putStrLn $ show gapArrow
--  putStrLn $ "perintahDOS Removed == " ++ unwords [
--    "genPDOSvert.hs ", atomOs, fromJust gapnya, spinnya, daftarFolder]
  generatedFATBAND <- genPBAND useOldBw spinnya "0" atomOs [daftarFolder]
  putStrLn $ show generatedFATBAND
  (xrangeatas,ticksbaru) <- genBandTicks foldernya
  putStrLn $ "========" ++ show xrangeatas ++ "=====" ++ show ticksbaru
  putStrLn $ "=== FINISHED PROCESSING : " ++ daftarLengkap
  plotSinglePic sisa useOldBw atomOs (colorId+1) ((xrangeatas,ticksbaru,arrow,generatedPBAND,generatedFATBAND):res)

genArrow :: [Char] -> [Char] -> [Char] -> [Char] -> Maybe String
genArrow _ _ _ _ = Nothing
genArrow' :: [Char] -> [Char] -> [Char] -> [Char] -> Maybe String
genArrow' bandGap valBandX valBandY condBandY =
        if (bandGap == "0")
           then Nothing
           else --(,) (Just $ "0:" ++ bandGap ++ "@" ++ valBandY)
                    (Just $ unlines [ concat ["set label sprintf ('{/Symbol D}=%.2feV',"
                                      , bandGap
                                      , ") at (" ++ valBandX ++ "-0.55),"
                                      ,"(" ++ valBandY ++ "-0.4) font ',12'"
                                      ]
                             , concat [ "set arrow from "
                                      , valBandX ++ "," ++ condBandY
                                      , " to " ++ valBandX ++ "," ++ valBandY
                                      , " heads noborder lw 2 lc rgb 'blue'"
                                      ]
                             ]
                    )

genPBAND :: Bool -> String -> String -> String -> [String] -> IO T.Text
--genPBAND oldBw spin invStat atomNos foldernya = do
genPBAND   _      _   _       []      _         = return ""
genPBAND   useOldBw  _   invStat atomNos foldernya = do
    let daftaratomOs =  map (splitOn "@") $ splitOn "-" atomNos
        daftarJudulSpinFolders = filter (/=[""]) $ map (splitOn "@") $ splitOn ":" $ unwords foldernya
        daftarAOJSF = zip ([1..] :: [Integer]) $ concat $ map (\a -> zip daftaratomOs $ repeat a) daftarJudulSpinFolders

    let err'' = map ( \(_,([o,a],[_,s,folder'])) -> concat [ "cd ", folder', "; BandWeight.py PROCAR.",cekSpin s invStat "1" "UP" "DN" , " ",a, " ", o]) daftarAOJSF
    if (useOldBw == True) then mapM_ system_ err''
                      else return ()

{--
  (\$2>6.08991818?0.4+\$2:\$2):(\$2>6.08991818&&\$2<6.09991818?1/0:\$3)
--}
    return $ T.pack $ concat $ intercalate [", "]
      $ map ( \(i,([o,a],[j,s,folder])) -> ["'",folder,"/bw.",a,".",o,".PROCAR.",cekSpin s invStat "1" "UP" "DN" ,".dat' u ($1>6.09971818?0.4+$1:$1):2:($1>6.08991818&&$1<6.09991818?0:$3*3) ls ",show i," ps variable title '",concat $ intersperse "." [o,a,j],"'"] :: [String] )
      $ daftarAOJSF

genPlotPlate :: Int -> T.Text -> T.Text
genPlotPlate colorId bnd = T.concat [ bnd , T.pack $ unwords [" u ($2>6.08991818?0.4+$2:$2):($2>6.08991818&&$2<6.09991818?1/0:$3) with line lc", show colorId, "title ''"]]

runGAPband :: [T.Text] -> String
runGAPband bndSpin =
  let (valencePoint, conductionPoint) =
        (\(a,b) -> (last a,head b)) $
        break (\[_,a] -> a > 0) $
        Data.List.sortBy (\[_,a] [_,b] -> if a < b then LT else GT ) $
        toLists $
        (¿ [1,2]) . fromLists $
        filter (/=[]) $  -- kolom $2 x (k-point) -- $3 y Energy
        map ( (map fst) . rights . (map T.double) . T.words ) $
        filter (not . T.isPrefixOf "#") $
        bndSpin

  in if (last conductionPoint - last valencePoint < 0.01)  then "0 0 0"
                                                           else unwords $ map show  [head valencePoint, head conductionPoint, last conductionPoint - last valencePoint]

cekSpin :: (Eq a, IsString a) => a -> a -> a-> p -> p -> p
cekSpin s "0" ok sOk sNo = if s == ok then sOk else sNo
cekSpin s _ ok sOk sNo = if s == ok then sNo else sOk

genBandTicks :: String -> IO (String,String)
genBandTicks foldernya = do
    lTicks'' <- inshell2text $ unwords ["tail -n2  ", foldernya ++ "/bnd*spin1|sed -e '/^ $\\|==/d'|sort -u|awk '{print $2}' "]
    namaTicks'' <- inshell2text $ unwords ["cat ",foldernya ++ "/syml.* "]
    let lTicks = (:) "0.0" $ map T.unpack $ filter (not . T.null) lTicks''
    let namaTicks =  map (\a -> if a == "Gamma" then "{/Symbol G}" else a) $ (map head $ filter (/=[]) $ map (snd . (splitAt 7) . words . T.unpack)  namaTicks'')
    return $ (,) (last lTicks)
           $ unwords $ intersperse "," $ map (\(a,b) -> unwords ["'"++a++"'",b]) $ zip  namaTicks lTicks

genTEMPGLT :: String -> String -> String -> String -> String -> String
           -> [String]
           -> String
genTEMPGLT tempDir judulUtama yr xrangeatas ticksbaru arrow plotplate =
    let subTitles = splitOn "#" judulUtama
     in unlines [ "#!/home/aku/bin/gnuplot -persist"
                , "reset"
                , "set term post landscape enhanced color 'Times-Roman' 12"
                , "set output '" ++ tempDir ++ "/tmp.ps'"
                , "rydberg=13.605"
                , "if (!exists('MP_LEFT'))   MP_LEFT = .1"
                , "if (!exists('MP_RIGHT'))  MP_RIGHT = .95"
                , "if (!exists('MP_BOTTOM')) MP_BOTTOM = .1"
                , "if (!exists('MP_TOP'))    MP_TOP = .9"
                , "if (!exists('MP_GAP'))    MP_GAP = 0.05"
                , unwords [ "set multiplot layout 1," ++ (show $ length subTitles )
                          , "margins screen MP_LEFT, MP_RIGHT, MP_BOTTOM, MP_TOP"
                          , "spacing screen MP_GAP"
                          ]
                , "set xzeroaxis"
                , "set grid"
                , "set key right top Left"
                , "set key samplen 1 spacing 1"
                , "set size ratio 1.5"
                , "set mytics 10"
                , "set xlabel 'Wave Vector'"
                , "set ylabel 'Energy-E_F (eV)'"
--                , "set title '" ++ judulUtama ++ "'"
                , "set yrange ["++ yr ++ "]"
                , "set xrange [0.0:" ++ xrangeatas ++ "]"
                , "set grid noy"
                , "set grid xtics lt 0 lc rgb 'black'"
                , "set style line 1 lt 2 lw 1 lc rgb '#e41a1c'"
                , "set style line 2 lt 6 lw 1 lc rgb '#377eb8'"
                , "set style line 3 lt 6 lw 1 lc rgb '#4daf4a'"
                , "set style line 4 lt 6 lw 1 lc rgb '#984ea3'"
                , "set style line 5 lt 6 lw 1 lc rgb '#ff7f00'"
                , "set style line 8 lt 1 lw 2 lc rgb '#ffcc99'"
                , "set style line 9 lt 1 lw 2 lc rgb '#808080'"
                , "set style line 10 lt 1 lw 2 lc rgb '#94ffb5'"
                , "set style line 11 lt 1 lw 2 lc rgb '#8f7c0'"
                , "set style line 12 lt 1 lw 2 lc rgb '#9dcc0'"
                , "set style line 13 lt 1 lw 2 lc rgb '#c2088'"
                , "set style line 14 lt 1 lw 2 lc rgb '#03380'"
                , "set style line 15 lt 1 lw 2 lc rgb '#ffa45'"
                , "set style line 16 lt 1 lw 2 lc rgb '#ffa8bb'"
                , "set style line 17 lt 1 lw 2 lc rgb '#42660'"
                , "set style line 18 lt 1 lw 2 lc rgb '#ff010'"
                , "set style line 19 lt 1 lw 2 lc rgb '#5ef1f2'"
                , "set style line 20 lt 1 lw 2 lc rgb '#0998f'"
                , "set style line 21 lt 1 lw 2 lc rgb '#e0ff66'"
                , "set style line 22 lt 1 lw 2 lc rgb '#74aff'"
                , "set style line 23 lt 1 lw 2 lc rgb '#9900'"
                , "set style line 24 lt 1 lw 2 lc rgb '#ffff80'"
                , "set style line 25 lt 1 lw 2 lc rgb '#ffff0'"
                , "set style line 26 lt 1 lw 2 lc rgb '#ff505'"
                , "set style arrow 1 heads size screen 0.01,90 lw 2 lc rgb 'navy'"
                , "set key bottom left Left"
                , "set xtics (" ++ ticksbaru ++ ")"
                , (unlines
                     $ intersperse noMidLabel
                     $ zipWith (\s p -> unlines $ ("set title '" ++ s ++ "'"):
                                        unwords [arrow, plotInit ,p]:
                                        [])
                       subTitles plotplate
                  )
                , endMultiplot
                ]

noMidLabel :: String
noMidLabel = unlines $ "unset ylabel":
                       []

data Opts = Opts {
    _fOut       :: String,
    _useOldBw   :: Bool,
    _judulUtama :: String,
    _yr         :: String,
    _atomOs     :: String,
    _targetDir  :: [String]
                 } deriving Show

optsParser :: Parser Opts
optsParser = Opts
             <$> strOption (long "out" <> short 'o' <> metavar "OUTPUT" <>
                          help "target output file" <> value "test")
             <*> switch (long "oldbw" <> short 'b'
                            <> help "Using old BandWeight files")
             <*> strOption (long "titles" <> short 't' <>
                            metavar "SUBTITLES" <>
                            help "Titles of each Figures" <>
                            value "")
             <*> strOption  ( long "yrange" <> short 'y' <>
                            metavar "MIN:MAX" <>
                            help "Y-range of the figure, ex. \"-2.5:7.3\"")
             <*> strOption (long "orbitals" <> short 'r' <> metavar "ORBITALS"
                            <> help "List of atomic orbitals" <> value "")
             <*> (many $ argument str (metavar "TARGETDIRS"))

{-
         "~/kanazawa/dev/cpvoh/app/f2 -o outputfile -b -y \"-9.5:2\" ext.nico2o4.normal/nico2o4.6G10:1:\"11G20.Majority\"#ext.nico2o4.normal/nico2o4.6G10:2:\"11G20.Minority\" ext.nico2o4.normal/nico2o4.1G0:1:\"G0.Majority\"#ext.nico2o4.normal/nico2o4.1G0:2:\"G0.Minority\"":
-}
withHelp :: ParserInfo Opts
withHelp = info
  (helper <*> optsParser)
  (fullDesc <> (progDesc
    "Band Structure Generator for ecalj, please make sure that TARGETDIRS are the last arguments."
  ))

max_stars_repo_name: hasanalrasyid/cpvoh

max_issues_repo_name: hasanalrasyid/cpvoh

max_forks_repo_name: hasanalrasyid/cpvoh


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

module Statistics.LinearRegression (
    -- * Simple linear regression functions
    linearRegression,
    linearRegressionRSqr,
    linearRegressionTLS,
    -- * Related functions
    correl,
    covar,
    -- * Estimated errors and distribution parameters
    linearRegressionMSE,
    linearRegressionDistributions,
    -- * Robust linear regression
    robustFit,
    nonRandomRobustFit,
    robustFitRSqr,
    -- ** Related types
    EstimationParameters(..),
    ErrorFunction,
    Estimator,
    EstimatedRelation,
    -- ** Provided values
    defaultEstimationParameters,
    linearRegressionError,
    linearRegressionTLSError,
    -- ** Helper functions
    converge,
    -- * References
    -- $references
    ) where

import qualified Data.Vector.Generic as G
import qualified Data.Vector.Unboxed as U
import Data.Vector.Generic (Vector, (!))
import Safe (at)
import System.Random
import System.Random.Shuffle (shuffleM)
import Control.Monad.Random.Class
import Control.Monad.Random (evalRand)
import Control.Monad (liftM)
import Data.Function (on)
import Data.List (minimumBy, sortBy)
import Data.Maybe (fromMaybe)
import qualified Statistics.Sample as S
import qualified Statistics.Distribution as D
import qualified Statistics.Distribution.Transform as T
import qualified Statistics.Distribution.StudentT as ST

--- * Simple linear regression

-- | Covariance of two samples
covar :: Vector v Double => v Double -> v Double -> Double
covar xs ys = covar' m1 m2 n xs ys
    where
          !n = fromIntegral $ G.length xs
          !m1 = S.mean xs
          !m2 = S.mean ys
{-# INLINE covar #-}

-- internal function that avoids duplicate calculation of means and lengths where possible
-- Note: trying to make the calculation even more efficient by subtracting m1*m1*n instead of individual subtractions increased errors, probably due to rounding issues.
covar' :: Vector v Double => Double -> Double -> Double -> v Double -> v Double -> Double
covar' m1 m2 n xs ys = G.sum (G.zipWith (*) (G.map (subtract m1) xs) (G.map (subtract m2) ys)) / (n-1)
{-# INLINE covar' #-}

-- | Pearson's product-moment correlation coefficient
correl :: Vector v Double => v Double -> v Double -> Double
correl xs ys = let !c = covar xs ys
                   !sx = S.stdDev xs
                   !sy = S.stdDev ys
               in c / (sx * sy)
{-# INLINE correl #-}

-- | Simple linear regression between 2 samples.
--   Takes two vectors Y={yi} and X={xi} and returns
--   (alpha, beta, r*r) such that Y = alpha + beta*X
--   and where r is the Pearson product-moment correlation
--   coefficient
linearRegressionRSqr :: Vector v Double => v Double -> v Double -> (Double, Double, Double)
linearRegressionRSqr xs ys = (alpha, beta, r2)
    where
          !c                   = covar' m1 m2 n xs ys
          !r2                  = c*c / (v1*v2)
          !(m1,v1)             = S.meanVarianceUnb xs
          !(m2,v2)             = S.meanVarianceUnb ys
          !n                   = fromIntegral $ G.length xs
          !beta                = c / v1
          !alpha               = m2 - beta * m1
{-# INLINE linearRegressionRSqr #-}

-- | Simple linear regression between 2 samples.
--   Takes two vectors Y={yi} and X={xi} and returns
--   (alpha, beta) such that Y = alpha + beta*X
linearRegression :: Vector v Double => v Double -> v Double -> (Double, Double)
linearRegression xs ys = (alpha, beta)
    where
        (alpha, beta, _) = linearRegressionRSqr xs ys
{-# INLINE linearRegression #-}

-- | The error (or residual) mean square of a sample w.r.t. an estimated regression line.
--   This serves as an estimate for the variance of the sampled data.
--   Accepts the regression parameters (alpha,beta) and the sample vectors X and Y.
linearRegressionMSE :: (Vector v Double, Vector v (Double, Double)) => (Double,Double) -> v Double -> v Double -> Double
linearRegressionMSE ab xs ys = (G.sum . G.map (linearRegressionError ab) . G.zip xs $ ys)/(n-2)
    where
        !n = fromIntegral $ G.length xs

-- | The estimated distributions of the regression parameters (alpha and beta) assuming normal, identical distributions of Y, the sampled data.
-- These can serve to get confidence intervals for the regression parameters.
-- Accepts the regression parameters (alpha,beta) and the sample vectors X and Y.
-- The distributions are StudnetT distributions centered at the estimated (alpha,beta) respectively, with parameter numbers n-2 (where n is the initial sample size) and with standard deviations that are extracted from the sampled data based on its MSE. See chapter 2 of reference [3] for details.
linearRegressionDistributions :: (Vector v Double, Vector v (Double, Double)) => (Double,Double) -> v Double -> v Double -> (T.LinearTransform ST.StudentT,T.LinearTransform ST.StudentT)
linearRegressionDistributions (alpha,beta) xs ys = (ST.studentTUnstandardized (n-2) alpha va,ST.studentTUnstandardized (n-2) beta vb)
    where
        !n = fromIntegral $ G.length xs
        !mse = linearRegressionMSE (alpha,beta) xs ys
        !vb = mse/(xv)
        !mx = S.mean xs
        !va = mse*(1/n+mx^2/xv)
        !xv = G.sum . G.map (\x -> (x-mx)^2) $ xs

-- | Total Least Squares (TLS) linear regression.
-- Assumes x-axis values (and not just y-axis values) are random variables and that both variables have similar distributions.
-- interface is the same as 'linearRegression'.
linearRegressionTLS :: Vector v Double => v Double -> v Double -> (Double,Double)
linearRegressionTLS xs ys = (alpha, beta)
    where
          !c                   = covar' m1 m2 n xs ys
          !b                   = (v1 - v2) / c
          !(m1,v1)             = S.meanVarianceUnb xs
          !(m2,v2)             = S.meanVarianceUnb ys
          !n                   = fromIntegral $ G.length xs
          !betas               = [(-b - sqrt(b^2+4))/2,(-b + sqrt(b^2+4)) /2]
          !beta                = if c > 0 then maximum betas else minimum betas
          !alpha               = m2 - beta * m1
{-# INLINE linearRegressionTLS #-}

-- | An estimated linear relation between 2 samples is (alpha,beta) such that Y = alpha + beta*X.
type EstimatedRelation = (Double,Double)

-- | An 'Estimator' is a function that generates an estimated linear regression based on 2 samples. This module provides two estimator functions:
-- 'linearRegression' and 'linearRegressionTLS'
type Estimator = (S.Sample -> S.Sample -> EstimatedRelation)

-- | An 'ErrorFunction' is a function that computes the error of a given point from an estimate. This module provides two error functions correspoinding to the two 'Estimator' functions it defines:
--
-- * Vertical distance squared via 'linearRegressionError' that should be used with 'linearRegression'
--
-- * Total distance squared vie 'linearRegressionTLSError' that should be used with 'linearRegressionTLS'
type ErrorFunction = (EstimatedRelation -> (Double,Double) -> Double)

-- | The robust fit algorithm used has various parameters that can be specified using the 'EstimationParameters' record.
data EstimationParameters = EstimationParameters {
    -- | Maximal fraction of outliers expected in the sample (default 0.25)
    outlierFraction     :: !Double,
    -- | Number of concentration steps to take for initial evaluation of a solution (default 3)
    shortIterationSteps :: !Int,
    -- | Maximal number of sampled subsets (pairs of points) to use as starting points (default 500)
    maxSubsetsNum       :: !Int,
    -- | If the initial sample is large, and thus gets subdivided, this is the number of candidate-estimations to take from each subgroup, on which complete convergence will be executed (default 10)
    groupSubsets        :: !Int,
    -- | Maximal size of sample that can be analyzed without any sub-division (default 600)
    mediumSetSize       :: !Int,
    -- | Maximal size of sample that does not require two-step sub-division (see reference article) (default 1500)
    largeSetSize        :: !Int,
    -- | Estimator function to use (default linearRegression)
    estimator           :: Estimator,
    -- | ErrorFunction to use (default linearRegressionError)
    errorFunction       :: ErrorFunction
    }

-- | Default set of parameters to use (see reference for details).
defaultEstimationParameters = EstimationParameters {
    outlierFraction = 0.25,
    shortIterationSteps = 3,
    maxSubsetsNum = 500,
    groupSubsets = 10,
    mediumSetSize = 600,
    largeSetSize = 1500,
    estimator = linearRegression,
    errorFunction = linearRegressionError
}

-- | linearRegression error function is the square of the /vertical/ distance of a point from the line.
linearRegressionError :: ErrorFunction
linearRegressionError (alpha,beta) (x,y) = (y-(beta*x+alpha))^2

-- | linearRegressionTLS error function is the square of the /total/ distance of a point from the line.
linearRegressionTLSError :: ErrorFunction
linearRegressionTLSError (alpha,beta) (x,y) = ey/(1+beta^2)
    where
        ey = linearRegressionError (alpha,beta) (x,y)

-- | Helper function to calculate the minimal expected size of uncontaminated data based on the maximal fraction of outliers.
setSize :: Vector v Double => EstimationParameters -> v Double -> Int
setSize ep xs = max (n `div` 2 + 1) . round $ (1-outlierFraction ep) * (fromIntegral n)
    where
        n = G.length xs

-- | Helper function that, given an initial estimated relation and the error of the perivous estimation, performs a "concentration" step, generating a new estimate based on a fraction of points laying closest to the previous estimate and estimates the error of the previous estimate based on the same fraction.
-- The result is an estimate that is at least as good as the previous one.
-- The reason the error is calculated for the previous parameters is calculation optimization.
concentrationStep :: Vector v Double => EstimationParameters -> v Double -> v Double -> (EstimatedRelation, Double) -> (EstimatedRelation, Double)
concentrationStep ep xs ys (prev, prev_err) = (new_estimate, new_err)
    where
        set_size = setSize ep xs
        xyerrors = map (\p -> (p,errorFunction ep prev p)) $ zip (G.toList xs) (G.toList ys)
        (xys,errors) = unzip . take set_size . sortBy (compare `on` snd) $ xyerrors
        (good_xs,good_ys) = unzip xys
        new_estimate = estimator ep (G.fromList good_xs) (G.fromList good_ys)
        new_err = sum errors

-- | Infinite set of consecutive concentration steps.
concentration :: Vector v Double => EstimationParameters -> v Double -> v Double -> EstimatedRelation -> [(EstimatedRelation, Double)]
concentration ep xs ys params = tail $ iterate (concentrationStep ep xs ys) (params,-1)

-- | Calculate the optimal (local minimum) estimate based on an initial estimate.
-- The local minimum may not be the global (a.k.a. best) estimate but starting from enough different initial estimates should yield the global optimum eventually.
converge :: Vector v Double => EstimationParameters -> v Double -> v Double -> EstimatedRelation -> EstimatedRelation
converge ep xs ys = fst . findConvergencePoint . concentration ep xs ys

-- | The convergence point is defined as the point the error estimate of which is equal to the next estimate's error.
findConvergencePoint :: Ord a => [(b,a)] -> (b,a)
findConvergencePoint (x:y:ys)
    | snd x <= snd y = x -- rounding issues my cause an actual increase in error resulting in an infinite loop so the actual stop condition is when the errors stop decreasing
    | otherwise = findConvergencePoint (y:ys)
findConvergencePoint xs = error "Too short a list for conversion (size < 2)"

-- | Many times there is no need for full concentration as bad initial estimates can be discovered after only a few concentration steps.
concentrateNSteps :: Vector v Double => EstimationParameters -> v Double -> v Double -> EstimatedRelation -> (EstimatedRelation,Double)
concentrateNSteps ep xs ys params = concentration ep xs ys params !! shortIterationSteps ep

-- | Finding a robust fit linear estimate between two samples. The procedure requires randomization and is based on the procedure described in the reference.
robustFit :: (MonadRandom m, Vector v Double) => EstimationParameters -> v Double -> v Double -> m EstimatedRelation
robustFit ep xs ys = do
    let n = G.length xs
-- For optimal performance the exact procedure executed depends on the set size.
    if n < 2
        then
            error "cannot fit an input of size < 2"
        else if n == 2
            then return $ lineParams ((G.head xs,G.head ys),(G.last xs,G.last ys))
            else
                liftM (candidatesToWinner ep xs ys) $ if n < mediumSetSize ep
                    then
                         singleGroupFitCandidates ep Nothing xs ys
                    else if n < largeSetSize ep
                        then largeGroupFitCandidates ep xs ys
                        else do
                            (nxs,nys) <- liftM unzip $ randomSubset (zip (G.toList xs) (G.toList ys)) (largeSetSize ep)
                            largeGroupFitCandidates ep (U.fromList nxs) (G.fromList nys)

-- | Robust fit yielding also the R-square value of the \"clean\" dataset.
robustFitRSqr :: (MonadRandom m, Vector v Double, Vector v (Double, Double)) => EstimationParameters -> v Double -> v Double -> m (EstimatedRelation,Double)
robustFitRSqr ep xs ys = do
    er <- robustFit ep xs ys
    let (good_xs,good_ys) = U.unzip . G.fromList . take (setSize ep xs) . sortBy (compare `on` errorFunction ep er) . G.toList $ G.zip xs ys
    return (er,correl good_xs good_ys ^ 2)

-- | A wrapper that executes 'robustFit' using a default random generator (meaning it is only pseudo-random)
nonRandomRobustFit :: Vector v Double => EstimationParameters -> v Double -> v Double -> EstimatedRelation
nonRandomRobustFit ep xs ys = evalRand (robustFit ep xs ys) (mkStdGen 1)

-- | Given a set of initial estimates converge them all and find the optimal one.
candidatesToWinner :: Vector v Double => EstimationParameters -> v Double -> v Double -> [EstimatedRelation] -> EstimatedRelation
candidatesToWinner ep xs ys = fst . minimumBy (compare `on` snd) . map (findConvergencePoint . concentration ep xs ys)

-- | for a large initial sample - subdivide it, then get candidates from each subgroup. Perform full convergence on all the candidates and return the best ones.
largeGroupFitCandidates :: (MonadRandom m, Vector v Double) => EstimationParameters -> v Double -> v Double -> m [EstimatedRelation]
largeGroupFitCandidates ep xs ys = do
    let n = G.length xs
    let sub_groups_num = n `div` (mediumSetSize ep `div` 2)
    let sub_groups_size = n `div` sub_groups_num
    shuffled <- shuffleM $ zip (G.toList xs) (G.toList ys)
    let sub_groups = map (G.unzip . U.fromList) $ splitTo sub_groups_size shuffled
    let sub_groups_candidates = maxSubsetsNum ep `div` sub_groups_num
    candidates_list <- mapM (applyTo $ singleGroupFitCandidates ep (Just sub_groups_candidates)) sub_groups
    let candidates = concat candidates_list
    return . map fst . take (groupSubsets ep) . sortBy (compare `on` snd) . map (findConvergencePoint . concentration ep xs ys) $ candidates

-- | For a single group (a group that will not be subdivided) pick an initial set of pairs of points, run a few steps on each, then return the most promising candidates.
singleGroupFitCandidates :: (MonadRandom m, Vector v Double) => EstimationParameters -> Maybe Int -> v Double -> v Double -> m [EstimatedRelation]
singleGroupFitCandidates ep m_subsets xs ys = do
    let all_pairs = allPairs $ zip (G.toList xs) (G.toList ys)
    let return_size = fromMaybe (maxSubsetsNum ep) m_subsets
    initial_sets <- if return_size > length all_pairs
        then return all_pairs
        else randomSubset all_pairs return_size
    return . map fst . take (groupSubsets ep) . sortBy (compare `on` snd) . map (concentrateNSteps ep xs ys . lineParams) $ initial_sets

-- | Find the line passing between two points. This is the initial estimate to use given two random points.
lineParams :: ((Double,Double),(Double,Double)) -> EstimatedRelation
lineParams ((x1,y1),(x2,y2)) = (alpha,beta)
    where
        beta = (y2-y1)/(x2-x1)
        alpha = y1 - beta*x1

-- | A list of all possible two-element pairs from a list.
allPairs :: [a] -> [(a,a)]
allPairs [] = []
allPairs [x] = []
allPairs [x,y] = [(x,y)]
allPairs (x:xs) = (zip xs . repeat $ x) ++ allPairs xs

-- | Get a random subset of a given size.
randomSubset :: MonadRandom m => [a] -> Int -> m [a]
randomSubset xs size = liftM (take size) $ shuffleM xs

-- | Split a list into sublists of length n.
splitTo :: Int -> [a] -> [[a]]
splitTo n = map (take n) . takeWhile (not . null) . iterate (drop n)

-- | Helper function to adjust parameter handling
applyTo :: (a->b->c) -> (a,b) -> c
applyTo f (x,y) = f x y

-- $references
--
-- * Two Dimensional Euclidean Regression (Stein) <http://www.dspcsp.com/pubs/euclreg.pdf>
--
-- * Computing LTS Regression For Large Data Sets (Rousseeuw and Driessen) <http://agoras.ua.ac.be/abstract/Comlts99.htm>
--
-- * Applied linear statistical models (Kutner et al.)

max_stars_repo_name: alpmestan/statistics-linreg

max_issues_repo_name: alpmestan/statistics-linreg

max_forks_repo_name: alpmestan/statistics-linreg


text:
module Stats (nintyFifth,  Estimate(..)) where

import qualified Data.Vector.Unboxed as U
import Statistics.Sample (mean)
import Statistics.Resampling (resample, fromResample)
import Statistics.Resampling.Bootstrap (bootstrapBCA, Estimate(..) )
import System.Random.MWC (create)

nintyFifth :: [Double] -> IO Estimate
nintyFifth sample = do
          g <- create
          resamples <- resample g [mean] 10000 sampleU -- (length sample^2) sampleU
--          print $ U.length $ fromResample $ head $ resamples
--          print resamples
          return $ head $ bootstrapBCA 0.95 sampleU [mean] resamples
  where
    sampleU = U.fromList sample

max_stars_repo_name: ku-fpg/ldpc

max_issues_repo_name: ku-fpg/ldpc

max_forks_repo_name: ku-fpg/ldpc


text:
module StateGenTests
  ( testStateGen
  , testStateArb
  , arbPropPure
  , arbPropDensityMatrix
  ) where

import Data.Complex (magnitude)
import HABQTlib.Data
import HABQTlib.RandomStates
import Numeric.LinearAlgebra (norm_2, sumElements, takeDiag, tr)
import qualified System.Random.MWC as MWC
import qualified Test.QuickCheck as QC
import Test.QuickCheck (Property)
import TestHelpers

class StateProp a where
  traceHermProp :: MWC.GenIO -> a -> Property

dimRankIO :: MWC.GenIO -> Dim -> IO (Dim, Rank)
dimRankIO gen ub = do
  dim <- MWC.uniformR (1, ub) gen
  r <- MWC.uniformR (1, dim) gen
  return (dim, r)

instance StateProp DensityMatrix where
  traceHermProp gen _ =
    QC.ioProperty $ do
      (dim, r) <- dimRankIO gen 100
      dm <- getDensityMatrix <$> genDM dim r
      return $
        (abs (1 - (magnitude . sumElements . takeDiag) dm) < 1e-12) &&
        (norm_2 (dm - tr dm) < 1e-12)

instance StateProp PureStateVector where
  traceHermProp gen _ =
    QC.ioProperty $ do
      (dim, _) <- dimRankIO gen 100
      sv <- getStateVector <$> genPureSV dim
      return $ abs (1 - norm_2 sv) < 1e-12

testStateGen :: IO ()
testStateGen = do
  gen <- MWC.createSystemRandom
  putStrLn "Testing density matrix generation:"
  QC.quickCheckWith
    QC.stdArgs {QC.maxSuccess = 10000}
    (traceHermProp gen :: DensityMatrix -> Property)
  putStrLn "Testing pure state generation:"
  QC.quickCheckWith
    QC.stdArgs {QC.maxSuccess = 10000}
    (traceHermProp gen :: PureStateVector -> Property)

arbPropPure :: PureStateVector -> Property
arbPropPure (PureStateVector sv) = QC.property $ abs (1 - norm_2 sv) < 1e-12

arbPropDensityMatrix :: DensityMatrix -> Property
arbPropDensityMatrix (DensityMatrix dm) =
  QC.property $
  (abs (1 - (magnitude . sumElements . takeDiag) dm) < 1e-12) &&
  (norm_2 (dm - tr dm) < 1e-12)

testStateArb :: IO ()
testStateArb = do
  putStrLn "Testing density matrix arbitrary:"
  QC.quickCheckWith
    QC.stdArgs {QC.maxSuccess = 1000}
    (arbPropDensityMatrix :: DensityMatrix -> Property)
  putStrLn "Testing pure state arbitrary:"
  QC.quickCheckWith
    QC.stdArgs {QC.maxSuccess = 10000}
    (arbPropPure :: PureStateVector -> Property)

max_stars_repo_name: Belinsky-L-V/HABQT

max_issues_repo_name: Belinsky-L-V/HABQT

max_forks_repo_name: Belinsky-L-V/HABQT


text:
module Main where

import System.IO
import Graphics.Gnuplot.Simple
import Numeric.LinearAlgebra

data DefIntegral a = DefIntegral (a->a) a a

--instance Show DefIntegral where
--  show (DefIntegral _ start end) = "Integral from" ++ show a ++ "to" ++ show b

--mixed :: (Enum a, Field a) => DefIntegral a -> a -> a -> a
--mixed (DefIntegral f start end) alpha steps =
--  sum $ map (\x -> x) (intervals start end steps)

aitken :: (Enum a, Field a, RealFrac a, Num a)
  => DefIntegral a
  -> a
  -> (DefIntegral a -> a -> a -> a)
  -> a
  -> (a, (a, a, a))
aitken integral@(DefIntegral _ a b) alpha method initial_step =
  let [s1, s2, s3] = map (method integral alpha) [1, 2, 4]
      error = 10**(-6)
      l = 2
      m = - log (abs $ (s3 - s2)/(s2 - s1)) / log l
      optimal_step = (b-a) * (abs (error * (1 - (b-a)**(-m) / (s2 - s1))))**(1/m)
  in (optimal_step, residue integral alpha method optimal_step)

residue1 :: (Enum a, Field a, RealFrac a, Ord a)
  => DefIntegral a
  -> a
  -> (DefIntegral a -> a -> a -> a)
  -> (a, a, a)
residue1 defintegral alpha method = residue defintegral alpha method 1

-- richardson method with r=2
-- returns tuple of integral value, residue and required step
residue :: (Enum a, Field a, RealFrac a, Ord a)
  => DefIntegral a
  -> a
  -> (DefIntegral a -> a -> a -> a)
  -> a
  -> (a, a, a)
residue (DefIntegral f start end) alpha method initial_step =
  let steps_from_step step = fromIntegral . floor $ (end - start) / step + 1
      quadrature step = method (DefIntegral f start end) alpha (steps_from_step step)
      find_m [s1, s2, s3] = - log (abs $ (s3 - s2)/(s2 - s1)) / log l
      error = 10**(-6)
      l = 2.0
      steps_ = [initial_step / l**i | i<-[0..]]

      iter steps =
        -- find vector (J, Cm, Cm+1, ...) and take J from it
        let steps' = (3><1) $ take 3 steps
            rhs = cmap quadrature steps'
            m = find_m . toList . flatten $ rhs
            lhs = (3><1) [1, 1, 1] ||| cmap (**m) steps' ||| cmap (**(m+1)) steps'
            integral = (lhs <\> rhs) `atIndex` (0, 0)
            res = integral - (rhs `atIndex` (2, 0))
        in (integral, res)
      find_residue steps =
        if abs res < error
        then (integral, res, steps !! 2)
        else find_residue (tail steps)
        where (integral, res) = iter steps
  in find_residue steps_

cardano :: (Field a) => a -> a -> a -> a -> [a]
cardano a b c d =
  let p = c/a - b**2/(3*a**2)
      q = 2*b**3 / (27*a**3) - b*c/(3*a**2) + d/a
      phi = acos(3*q/(2*p) * sqrt(-3/p))
      y1 = 2*sqrt(-p/3) * cos(phi/3)
      y2 = 2*sqrt(-p/3) * cos(phi/3 + 2*pi/3)
      y3 = 2*sqrt(-p/3) * cos(phi/3 - 2*pi/3)
  in [y1 - b/(3*a), y2 - b/(3*a), y3 - b/(3*a)]

gauss_helper :: (Enum a, Field a) => DefIntegral a -> a -> (a, a) -> a
gauss_helper (DefIntegral f start end) alpha (left, right) =
  let -- 0) substition
      left' = end - left
      right' = end - right
      avg' = (right' + left') / 2
      f' x = f $ end - x

      -- 1) find weights
      weight n = - (right'**(n-alpha+1) - left'**(n-alpha+1)) / (n-alpha+1)

      --2) solve linear system for a's
      ws_lhs = (3><3)
        [ weight 0, weight 1, weight 2
        , weight 1, weight 2, weight 3
        , weight 2, weight 3, weight 4]

      ws_rhs = (3><1)
        [ -weight 3
        , -weight 4
        , -weight 5]

      coefs_a = ws_lhs <\> ws_rhs

      -- 3) solve 3-order equation for xs
      [d, c, b] = toList . flatten $ coefs_a
      xs = (1><3) $ cardano 1 b c d

      -- 4) solve linear system for A's
      helper_lhs = (1><3) [1, 1, 1] === xs === cmap (**2) xs
      helper_rhs = (3><1)
        [ weight 0
        , weight 1
        , weight 2]
      coefs_A = helper_lhs <\> helper_rhs

      -- then sum Ai*fi
      fs = cmap f' xs
      result_vector = fs <> coefs_A
      [[result]] = toLists result_vector
  in result

gauss :: (Enum a, Field a) => DefIntegral a -> a -> a -> a
gauss integral@(DefIntegral f start end) alpha steps =
  sum $ map (gauss_helper integral alpha) (intervals start end steps)

general_nc_helper :: (Enum a, Field a) => Int -> DefIntegral a -> a -> (a, a) -> a
general_nc_helper order (DefIntegral f start end) alpha (left, right) =
  let  -- substitution
      left' = end - left
      right' = end - right
      f' x = f $ end - x

      xs = linspace order (end-left, end-right)
      fs = asColumn $ cmap f' xs

      weight n = - (right'**(n-alpha+1) - left'**(n-alpha+1)) / (n-alpha+1)

      weights = cmap weight $ (order><1) [0..toEnum order-1]

      xs_matrix = (order><order) [x^p | p<-[0..order-1], x<-toList xs]

      a_coefs = tr $ xs_matrix <\> weights
      result_vector = a_coefs <> fs

      [[result]] = toLists result_vector
  in result

newton_cotes_helper :: (Enum a, Field a) => DefIntegral a -> a -> (a, a) -> a
newton_cotes_helper (DefIntegral f start end) alpha (left, right) =
  let -- substitution
      left' = end - left
      right' = end - right
      avg' = (right' + left') / 2
      f' x = f $ end - x

      weight n = - (right'**(n-alpha+1) - left'**(n-alpha+1)) / (n-alpha+1)

      weights = (3><1) [weight 0, weight 1, weight 2]

      xs = (3><3)
        [ 1, 1, 1
        , left', avg', right'
        , left'**2, avg'**2, right'**2]

      fs = (3><1) [f' left', f' avg', f' right']

      a_coefs = tr $ xs <\> weights

      result_vector = a_coefs <> fs
  in result_vector `atIndex` (0, 0)

--alpha is power of p(x)
newton_cotes :: (Enum a, Field a) => DefIntegral a -> a -> a -> a
newton_cotes integral@(DefIntegral f start end) alpha steps =
  sum $ map (general_nc_helper 3 integral alpha) (intervals start end steps)

-- mixed with same degree=5
mixed1 :: (Enum a, Field a) => DefIntegral a -> a -> a -> a
mixed1 integral@(DefIntegral f start end) alpha steps =
  sum . map (\(interval, i) -> helper interval i) $ zip (intervals start end steps) [1..]
  where helper interval i = if odd i
                            then general_nc_helper 6 integral alpha interval
                            else gauss_helper integral alpha interval

-- mixed with different degrees: gauss with 5 and nc with 2
mixed2 :: (Enum a, Field a) => DefIntegral a -> a -> a -> a
mixed2 integral@(DefIntegral f start end) alpha steps =
  sum . map (\(interval, i) -> helper interval i) $ zip (intervals start end steps) [1..]
  where helper interval i = if odd i
                            then general_nc_helper 3 integral alpha interval
                            else gauss_helper integral alpha interval

left_sum :: (Fractional a, Num a, Enum a) => DefIntegral a -> a -> a
left_sum (DefIntegral f start end) steps =
  sum . map (\(left, right) -> (right - left) * (f left)) $ intervals start end steps

average_sum :: (Fractional a, Num a, Enum a) => DefIntegral a -> a -> a
average_sum (DefIntegral f start end) steps =
  sum . map (\(left, right) -> (right - left) * (f $ (left + right) / 2)) $ intervals start end steps

trapezoid_sum :: (Fractional a, Num a, Enum a) => DefIntegral a -> a -> a
trapezoid_sum (DefIntegral f start end) steps =
  sum . map (\(left, right) -> (right - left) * (f right + f left) / 2) $ intervals start end steps

simpson_sum :: (Fractional a, Num a, Enum a) => DefIntegral a -> a -> a
simpson_sum (DefIntegral f start end) steps =
  halfstep / 3 * (sum . map (\(left, middle, right) -> f left + 4 * f middle +  f right) $ interv)
  where step = (end - start) / steps
        halfstep = step / 2
        segment = [start + halfstep, start + 3*halfstep .. end - halfstep]
        interv = map (\x -> (x - halfstep, x, x + halfstep)) segment

simpsons :: (Fractional a, Num a) => DefIntegral a -> a
simpsons (DefIntegral f start end) =
  (end - start) / 6 * (f start + 4 * f ((start + end) / 2) + f end)

intervals :: (Fractional a, Num a, Enum a) => a -> a -> a -> [(a, a)]
intervals start end n_intervals =
  zip intervals' (tail intervals')
  where intervals' = [start, start + step .. end - step] ++ [end]
        -- add end to the end to escape problems with floating point arithmetics
        step  = (end - start) / n_intervals

tests :: IO ()
tests = do
  putStrLn "hello world"
  -- let a = 0 :: Double
  -- let b = 10
  -- let func = \x -> x
  let func = \x -> 3.7 * cos(1.5 * x) * exp(-4*x / 3) + 2.4 * sin(4.5 * x) * exp(2*x / 3) + 4
  let a = 1.8 :: Double
  let b = 2.3

  let integral = DefIntegral func a b
  let test_input = [2, 3, 5, 10, 50, 100, 200]

  let test method = map (method integral) test_input
  let test2 method = map (method integral 0.6) test_input

  putStrLn $ "steps: " ++ show test_input
  putStrLn $ "newton-cotes: " ++ show (test2 newton_cotes)
  putStrLn $ "gauss: " ++ show (test2 gauss)
  putStrLn $ "mixed1: " ++ show (test2 mixed1)
  putStrLn $ "mixed2: " ++ show (test2 mixed2)
  putStrLn "------------------------------\n"
  -- should be 1.18515
  putStrLn $ "left_sum: " ++ show (test left_sum)
  putStrLn $ "avg_sum: " ++ show  (test average_sum)
  putStrLn $ "trap_sum: " ++ show (test trapezoid_sum)
  putStrLn $ "simpsons: " ++ show (test simpson_sum)

  putStrLn $ "\n##########################\n"

  let residue_test method = residue1 integral 0.6 method

  putStrLn $ "residue newton_cotes: " ++ show (residue_test newton_cotes)
  putStrLn $ "residue gauss: " ++ show (residue_test gauss)
  putStrLn $ "residue mixed1: " ++ show (residue_test mixed1)
  putStrLn $ "residue mixed2: " ++ show (residue_test mixed2)

  let aitken_gauss_test = aitken integral 0.6 gauss 0.1
  let (opt_step, aitken_result) = aitken_gauss_test
  putStrLn $ "optimal step = " ++ show opt_step
  putStrLn $ "aitken gauss: " ++ show aitken_result

plots :: IO()
plots = do
  let func = \x -> 3.7 * cos(1.5 * x) * exp(-4*x / 3) + 2.4 * sin(4.5 * x) * exp(2*x / 3) + 4
  let a = 1.8 :: Double
  let b = 2.3
  let integral = DefIntegral func a b
  let steps = [3.0,4.0..500]

  let integral_value = 2.37880192961486
  --let integral_value = 1.18515

  let map_left = map (left_sum integral) steps
  let map_avg = map (average_sum integral) steps
  let map_trap = map (trapezoid_sum integral) steps
  let map_simps = map (simpson_sum integral) steps

  let map_nc = map (newton_cotes integral 0.6) steps
  let map_gauss = map (gauss integral 0.6) steps

  let left = zip steps $ map (-integral_value+) map_left
  let avg = zip steps $ map (-integral_value+) map_avg
  let trap = zip steps $ map (-integral_value+) map_trap
  let simps = zip steps $ map (-integral_value+) map_simps
  --plotList [] $ zip steps $ map (-integral_value+) map_nc
  --plotList [] $ zip steps $ map (-integral_value+) map_gauss

  let helper values name = (defaultStyle {lineSpec = CustomStyle [LineTitle name]}, values)
  plotListsStyle [Title "Comparacent of different quadrature methods"]
    [ helper left "left"
    , helper avg "avg"
    , helper trap "trapezoid"
    , helper simps "simpson"]

outputToFile :: IO()
outputToFile = do
  let func = \x -> 3.7 * cos(1.5 * x) * exp(-4*x / 3) + 2.4 * sin(4.5 * x) * exp(2*x / 3) + 4
  let a = 1.8 :: Double
  let b = 2.3
  let integral = DefIntegral func a b
  let steps = [3.0,4.0..500]

  let integral_value = 2.37880192961486

  let map_left = map (left_sum integral) steps
  let map_avg = map (average_sum integral) steps
  tests

  file <- openFile "outputnc.txt" WriteMode
  --hPutStrLn file ("integral_result = " ++ show integral_value ++ "\n")
  --hPutStrLn file ("steps = " ++ (show steps) ++ "\n")
  --hPutStrLn file ("left_sum = " ++ (show $ map (left_sum integral) steps) ++ "\n")
  --hPutStrLn file ("avg_sum = " ++ (show $ map (average_sum integral) steps) ++ "\n")
  --hPutStrLn file ("trap_sum = " ++ (show $ map (trapezoid_sum integral) steps) ++ "\n")
  hPutStrLn file ("simpson_sum = " ++ (show $ map (simpson_sum integral) steps) ++ "\n")
  --hPutStrLn file ("newton_cotes_sum = " ++ (show $ map (newton_cotes integral 0.6) steps) ++ "\n")
  hClose file

main :: IO ()
main = do
  --outputToFile
  tests

max_stars_repo_name: zykt/quadratures

max_issues_repo_name: zykt/quadratures

max_forks_repo_name: zykt/quadratures


text:
module Algebra where

import qualified Data.Text as T
import qualified Data.Map as M
import qualified Data.Vector as V
import Data.Maybe
import Statistics.LinearRegression

import Types

selectKeyData :: T.Text -> [Sample ItemsMap] -> [Sample Int]
selectKeyData key samples = [sample {sampleItems = select sample} | sample <- samples]
  where
    select sample = fromMaybe 0 $ M.lookup key (sampleItems sample)

prepareSystem :: [Sample Int] -> (V.Vector Double, V.Vector Double)
prepareSystem samples =
      (V.fromList $ map sampleTime samples,
       V.fromList $ map (fromIntegral . sampleItems) samples)

growCoefficient :: T.Text -> [Sample ItemsMap] -> Double
growCoefficient key samples =
  let (xs, ys) = prepareSystem $ selectKeyData key samples
      (_, beta) = linearRegression xs ys
  in  beta


text:
-- | Simple linear regression example for the README.

import Lib(someFunc)

import Control.Monad (replicateM, replicateM_, zipWithM)
import System.Random (randomIO)
import Test.HUnit (assertBool)

import qualified TensorFlow.Core as TF
import qualified TensorFlow.GenOps.Core as TF
import qualified TensorFlow.Gradient as TF
import qualified TensorFlow.Ops as TF

import qualified Data.Random.Normal as Normal
import qualified Statistics.Sample as Stat
import qualified Data.Vector as Vector
import Text.Printf

-- | R² of the model Ŷ = aX+b
-- Note that 'TF.gradient' doesn't support Double
linearModelRSquared :: [Float] -> [Float] -> Float -> Float
                    -> Float
linearModelRSquared xs ys a b = 1 - ssRes / ssTot
        where
                es    = zipWith (-) ys yHats
                yHats = [ a*x+b | x <- xs ]
                -- realToFrac is only ok in this situation
                yMean = realToFrac $ Stat.mean $ Vector.fromList $ (map realToFrac ys)
                ssTot = sum [ (y - yMean) ^ 2 | y <- ys ]
                ssReg = sum [ (y - yMean) ^ 2 | y <- yHats ]
                ssRes = sum [ e ^ 2 | e <- es ]

main :: IO ()
main = do
        -- Generate a list of Gaussian errors
        noise <- replicateM 100 (Normal.normalIO' (0, 0.1))

        let xData = tail [0,0.01..1]
            yData = [ 3*x + 8 + e | (x,e) <- zip xData noise]

       -- Fit the linear regression model.
        (w, b) <- fit xData yData

       -- Calculate and print the R-squared for the model
        let rSquared = linearModelRSquared xData yData w b
        printf "R^2 = %f\n" rSquared

fit :: [Float] -> [Float] -> IO (Float, Float)
fit xData yData = TF.runSession $ do
    -- Create tensorflow constants for x and y.
    let x = TF.vector xData
        y = TF.vector yData
    -- Create scalar variables for slope and intercept.
    w <- TF.initializedVariable 0
    b <- TF.initializedVariable 0
    -- Define the loss function.
    let yHat = (x `TF.mul` w) `TF.add` b
        loss = TF.square (yHat `TF.sub` y)
    -- Optimize with gradient descent.
    trainStep <- gradientDescent 0.001 loss [w, b]
    replicateM_ 1000 (TF.run trainStep)
    -- Return the learned parameters.
    (TF.Scalar w', TF.Scalar b') <- TF.run (w, b)
    return (w', b')

gradientDescent :: Float
                -> TF.Tensor TF.Build Float
                -> [TF.Tensor TF.Ref Float]
                -> TF.Session TF.ControlNode
gradientDescent alpha loss params = do
    let applyGrad param grad =
            TF.assign param (param `TF.sub` (TF.scalar alpha `TF.mul` grad))
    TF.group =<< zipWithM applyGrad params =<< TF.gradients loss params


text:
module Render(render) where

import Core
import Codec.Picture
import Codec.Picture.Types
import Control.Monad
import System.Random
import Control.Monad.Random

import Numeric.LinearAlgebra

deltas :: RenderSettings -> (Vec3, Vec3)
deltas settings  =
    let
        scale_x = 1 / (fromIntegral (width settings) - 1)
        scale_y = 1 / (fromIntegral (height settings) - 1)
        dx = scale scale_x (topRight settings - topLeft settings)
        dy = scale scale_y (bottomRight settings - topRight settings)
    in (dx,dy)

calc_ray :: RenderSettings -> Int -> Int -> Ray
calc_ray settings x y =
    let
        (dx,dy) = deltas settings
        point = topLeft settings + scale (fromIntegral x) dx + scale (fromIntegral y) dy
        dir = point - origin settings
    in (Ray point dir)

render :: RenderSettings -> Scene -> IO DynamicImage
render settings scene = do
    mut_img <- createMutableImage (width settings) (height settings) (background settings)
    let pixels = [(x,y) | x <- [1..width settings], y <- [1..height settings]]
    let (dx,dy) = deltas settings
    forM_ pixels $ \(x,y) -> do
        let (Ray point dir) = calc_ray settings x y
        colors <- forM [1..antialiasing settings] $ \_ -> do
            kx <- randomRIO (-1,1)
            ky <- randomRIO (-1,1)
            gen <- getStdGen
            let ray' = Ray (point) (dir + scale kx dx + scale ky dy)
            let (color, gen') = runRand (traceRay (background settings) scene ray') gen
            setStdGen gen'
            return $! (1/fromIntegral (antialiasing settings),color)
        writePixel mut_img (x-1) (y-1) (weigh colors)
    img <- unsafeFreezeImage mut_img
    return $! ImageRGBF img

max_stars_repo_name: danrocag/raytracer

max_issues_repo_name: danrocag/raytracer

max_forks_repo_name: danrocag/raytracer


text:
{-# LANGUAGE CPP #-}
#ifndef MIN_VERSION_base
#include "../dist/build/autogen/cabal_macros.h"
#endif
module Calico.Base (

    -- * Basic types

    -- ** Maybe
    Data.Maybe.Maybe(..)
  , Data.Maybe.maybe
  , Data.Maybe.isJust
  , Data.Maybe.isNothing
  , Data.Maybe.fromMaybe
  , Data.Maybe.listToMaybe
  , Data.Maybe.maybeToList
  , Data.Maybe.catMaybes
  , Data.Maybe.mapMaybe

    -- ** Either
  , Data.Either.Either(..)
  , Data.Either.either

    -- ** Bool
  , Data.Bool.Bool(..)
  , (Data.Bool.&&)
  , (Data.Bool.||)
  , Data.Bool.not
  , Data.Bool.otherwise

    -- ** Numeric types
#if MIN_VERSION_base(4, 8, 0)
  , Numeric.Natural.Natural
#endif
  , Prelude.Integer
  , Prelude.Int
  , Prelude.Rational
  , Prelude.Float
  , Prelude.Double

    -- ** Char
  , Data.Char.Char
  , Data.Char.toUpper
  , Data.Char.toLower
  , Data.Char.ord
  , Data.Char.chr

    -- ** Word
  , Data.Word.Word
  , Data.Word.Word8
  , Data.Word.Word16
  , Data.Word.Word32
  , Data.Word.Word64

    -- ** Ratio
  , Data.Ratio.Ratio
  , (Data.Ratio.%)
  , Data.Ratio.numerator
  , Data.Ratio.denominator
  , Data.Ratio.approxRational

    -- ** Complex
  , Data.Complex.Complex(..)
  , Data.Complex.cis
  , Data.Complex.conjugate
  , Data.Complex.imagPart
  , Data.Complex.magnitude
  , Data.Complex.mkPolar
  , Data.Complex.polar
  , Data.Complex.phase
  , Data.Complex.realPart

    -- ** List
  , (Data.List.!!)
  , Data.List.intercalate
  , Data.List.permutations
  , Data.List.zip
  , Data.List.zipWith
#if MIN_VERSION_base(4, 8, 0)
  , Data.List.uncons
#else
  , uncons
#endif
  , Data.List.unzip

    -- ** String
  , Data.String.String
  , Data.String.IsString
  , Text.Printf.printf
  , Text.Printf.PrintfType

    -- * Basic type classes

    -- ** Ordering
  , Data.Eq.Eq(..)
  , Data.Ord.Ordering(..)
  , Data.Ord.Ord(..)

    -- ** Numerics
  , Prelude.Enum(..)
  , Prelude.Bounded(..)
  , Prelude.Num(..)
  , Prelude.Real(..)
  , Prelude.Integral(..)
  , Prelude.Fractional(..)
  , Prelude.Floating(..)
  , Prelude.RealFrac(..)
  , Prelude.RealFloat(..)
  , Prelude.subtract
  , Prelude.even
  , Prelude.odd
  , Prelude.gcd
  , Prelude.lcm
  , (Prelude.^)
  , (Prelude.^^)
  , Prelude.fromIntegral
  , Prelude.realToFrac

    -- ** Bits
  , Data.Bits.Bits
  , (Data.Bits..&.)
  , (Data.Bits..|.)
  , Data.Bits.xor
  , Data.Bits.complement
  , Data.Bits.shift
  , Data.Bits.rotate
  , Data.Bits.bit
  , Data.Bits.setBit
  , Data.Bits.clearBit
  , Data.Bits.complementBit
  , Data.Bits.testBit
  , Data.Bits.isSigned
  , Data.Bits.shiftL
  , Data.Bits.shiftR
  , Data.Bits.rotateL
  , Data.Bits.rotateR

    -- ** Read
  , Text.Read.ReadS
  , Text.Read.Read
  , Text.Read.readsPrec
  , Text.Read.readList
  , Text.Read.reads
  , Text.Read.readParen
  , Text.Read.read
  , maybeRead
  , Text.Read.lex

    -- ** Show
  , Text.Show.ShowS
  , Text.Show.Show(..)
  , Text.Show.shows
  , Text.Show.showChar
  , Text.Show.showString
  , Text.Show.showParen
  , Text.Show.showListWith

    -- ** Monoid
  , Data.Monoid.Monoid(..)
  , (Data.Monoid.<>)

    -- ** Foldable
  , Data.Foldable.Foldable(..)
#if !MIN_VERSION_base(4, 8, 0)
  , Data.Foldable.toList
  , null
  , length
  , Data.Foldable.elem
  , Data.Foldable.maximum
  , Data.Foldable.minimum
  , Data.Foldable.sum
  , Data.Foldable.product
#endif
  , Data.Foldable.foldlM
  , Data.Foldable.foldrM
  , Data.Foldable.traverse_
  , Data.Foldable.for_
  , Data.Foldable.sequenceA_
  , Data.Foldable.asum
  , Data.Foldable.concatMap
  , Data.Foldable.and
  , Data.Foldable.or
  , Data.Foldable.any
  , Data.Foldable.all
  , sum'
  , prod'
  , Data.Foldable.maximumBy
  , Data.Foldable.minimumBy
  , Data.Foldable.notElem
  , Data.Foldable.find

    -- ** Traversable
  , Data.Traversable.Traversable(..)
  , Data.Traversable.for
  , Data.Traversable.forM
  , Data.Traversable.mapAccumL
  , Data.Traversable.mapAccumR
  , Data.Traversable.fmapDefault
  , Data.Traversable.foldMapDefault

    -- ** Category
  , Control.Category.Category
  , Control.Category.id
  , (Control.Category..)
  , (Control.Category.<<<)
  , (Control.Category.>>>)

    -- ** Arrow
  , Control.Arrow.Arrow
  , Control.Arrow.arr
  , Control.Arrow.first
  , Control.Arrow.second
  , (Control.Arrow.***)
  , (Control.Arrow.&&&)
  , Control.Arrow.returnA
  , (Control.Arrow.^>>)
  , (Control.Arrow.>>^)
  , Control.Arrow.ArrowChoice
  , Control.Arrow.left
  , Control.Arrow.right
  , (Control.Arrow.+++)
  , (Control.Arrow.|||)

    -- ** Functor
  , Data.Functor.Functor(..)
  , (Data.Functor.<$>)
#if MIN_VERSION_base(4, 7, 0)
  , (Data.Functor.$>)
#else
  , ($>)
#endif
  , Data.Functor.void

    -- ** Applicative
  , Control.Applicative.Applicative(..)
  , Control.Applicative.Alternative(..)
  , Control.Applicative.Const(..)
  , (Control.Applicative.<**>)
  , Control.Applicative.optional

    -- ** Monad
  , Control.Monad.Monad(..)
  , Control.Monad.MonadPlus(..)
  , (Control.Monad.=<<)
  , (Control.Monad.<=<)
  , (Control.Monad.>=>)
  , Control.Monad.join
  , Control.Monad.forever
  , Control.Monad.guard
  , Control.Monad.when
  , Control.Monad.unless

    -- ** Generic
  , GHC.Generics.Generic

    -- ** Typeable
  , Data.Typeable.Typeable

    -- ** Type-level programming
#if MIN_VERSION_base(4, 7, 0)
  , Data.Proxy.Proxy(Proxy)
  , Data.Proxy.asProxyTypeOf
  , Data.Proxy.KProxy(KProxy)
#endif

    -- * Functions

    -- ** Tuples
  , Data.Tuple.fst
  , Data.Tuple.snd
  , Data.Tuple.curry
  , Data.Tuple.uncurry
  , Data.Tuple.swap

    -- ** Utility
  , (Data.Function.$)
  , Data.Function.flip
  , Data.Function.fix
  , Prelude.seq
  , (Prelude.$!)
  , Prelude.asTypeOf
  , Data.Function.on
  , Prelude.until
  , Prelude.error

    -- ** Debugging
  , Prelude.undefined
  , Debug.Trace.trace
  , Debug.Trace.traceShow
  , Debug.Trace.traceStack
#if MIN_VERSION_base(4, 7, 0)
  , Debug.Trace.traceId
  , Debug.Trace.traceShowId
#else
  , traceId
  , traceShowId
#endif

  ) where
import qualified Prelude
import qualified Control.Applicative
import qualified Control.Arrow
import qualified Control.Category
import qualified Control.Monad
import qualified Data.Bits
import qualified Data.Bool
import qualified Data.Char
import qualified Data.Complex
import qualified Data.Either
import qualified Data.Eq
import qualified Data.Foldable
import qualified Data.Functor
import qualified Data.Function
import qualified Data.List
import qualified Data.Maybe
import qualified Data.Monoid
import qualified Data.Ord
#if MIN_VERSION_base(4, 7, 0)
import qualified Data.Proxy
#endif
import qualified Data.Ratio
import qualified Data.String
import qualified Data.Traversable
import qualified Data.Typeable
import qualified Data.Word
import qualified Data.Tuple
import qualified Debug.Trace
import qualified GHC.Generics
#if MIN_VERSION_base(4, 8, 0)
import qualified Numeric.Natural (Natural)
#endif
import qualified Text.Read
import qualified Text.Show
import qualified Text.Printf
import Data.Maybe (Maybe(..))
import Data.String (String)
import Text.Read (Read, reads)

#if !MIN_VERSION_base(4, 7, 0)
infixl 4 $>
($>) :: Functor f => f a -> b -> f b
($>) = Data.Functor.flip (Data.<$)
#endif

#if !MIN_VERSION_base(4, 8, 0)
-- | Note: if base is older than 4.8, this may not be very efficient.
null :: Data.Foldable.Foldable t => t a -> Data.Bool.Bool
null = Data.Foldable.foldr (\_ _ -> Data.Bool.False) Data.Bool.True

-- | Note: if base is older than 4.8, this will always be /O(n)/.  If
--   performance is critical, consider explicitly importing from
--   @Data.Foldable@ or use the specialized versions.
length :: Data.Foldable.Foldable t => t a -> Prelude.Int
length = Data.Foldable.foldr (\_ x -> Prelude.succ x) 0
#endif

-- | Strict version of 'sum'.
sum' :: (Data.Foldable.Foldable t, Prelude.Num a) => t a -> a
sum' = Data.Foldable.foldl' (Prelude.+) 0

-- | Strict version of 'prod'.
prod' :: (Data.Foldable.Foldable t, Prelude.Num a) => t a -> a
prod' = Data.Foldable.foldl' (Prelude.*) 1

#if !MIN_VERSION_base(4, 7, 0)
traceId :: Show a => a -> a
traceId a = Debug.Trace.trace a a

traceShowId :: Show a => a -> a
traceShowId a = Debug.Trace.traceShow a a
#endif

#if !MIN_VERSION_base(4, 8, 0)
uncons :: [a] -> Maybe (a, [a])
uncons []       = Nothing
uncons (x : xs) = Just (x, xs)
#endif

maybeRead :: Read a => String -> Maybe a
maybeRead s = do
  case reads s of
    [(x, "")] -> Just x
    _         -> Nothing

max_stars_repo_name: Rufflewind/calico-hs

max_issues_repo_name: Rufflewind/calico-hs

max_forks_repo_name: Rufflewind/calico-hs


text:
#!/usr/bin/env stack
-- stack runghc --resolver lts-16 --package ad --package hmatrix-vector-sized-0.1.3.0 --package vector-sized --package hmatrix --package comonad --package free

-- | Source file accompanying
-- https://blog.jle.im/entry/hamiltonian-dynamics-in-haskell.html

{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE RankNTypes          #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns        #-}

import           Data.Foldable
import           Data.Maybe
import           GHC.TypeLits
import           Numeric.AD
import           Numeric.LinearAlgebra.Static
import           Numeric.LinearAlgebra.Static.Vector
import qualified Control.Comonad                     as C
import qualified Control.Comonad.Cofree              as C
import qualified Data.Vector.Sized                   as V

-- | A @'System' m n@ represents a system parameterized by @n@ generalized
-- coordinates, moving in an "underlying" @m@-dimensional cartesian
-- coordinate system.
data System m n = System
    { sysInertia       :: R m                         -- ^ 'm' vector
    , sysCoords        :: R n -> R m                  -- ^ f
    , sysJacobian      :: R n -> L m n                -- ^ J_f
    , sysHessian       :: R n -> V.Vector n (L m n)   -- ^ H_f
    , sysPotential     :: R n -> Double               -- ^ U
    , sysPotentialGrad :: R n -> R n                  -- ^ grad U
    }

-- | A system's position in configuration space
data Config n = Config
    { confPositions  :: R n
    , confVelocities :: R n
    }
  deriving Show

-- | A system's position in phase space
data Phase n = Phase
    { phasePositions :: R n
    , phaseMomenta   :: R n
    }
  deriving Show

-- | Given a position in @n@-dimensional generalized coordinates, return
-- the position in the underlying @m@-dimensional cartesian coordinate
-- system.
underlyingPosition
    :: System m n
    -> R n
    -> R m
underlyingPosition = sysCoords

-- | Give the system's generalized momenum, given its position in
-- configuration space.
momenta
    :: (KnownNat n, KnownNat m)
    => System m n
    -> Config n
    -> R n
momenta s (Config q v) = tr j #> mHat #> j #> v
  where
    j    = sysJacobian s q
    mHat = diag (sysInertia s)

-- | Convert a position in configuration space to a position in phase space
toPhase
    :: (KnownNat n, KnownNat m)
    => System m n
    -> Config n
    -> Phase n
toPhase s c = Phase (confPositions c) (momenta s c)

-- | Shift around the Hessian into the shape expected
tr2 :: (KnownNat m, KnownNat n)
    => V.Vector m (L n n)
    -> V.Vector n (L m n)
tr2 = fmap rowsL . traverse lRows
{-# INLINE tr2 #-}

-- | Vector of vectors to matrix
vec2l
    :: KnownNat n
    => V.Vector m (V.Vector n Double)
    -> L m n
vec2l = rowsL . fmap gvecR
{-# INLINE vec2l #-}

-- | Make a system given its inertias, coordinate functions, and potential
-- energy function
mkSystem
    :: (KnownNat m, KnownNat n)
    => R m
    -> (forall a. RealFloat a => V.Vector n a -> V.Vector m a)
    -> (forall a. RealFloat a => V.Vector n a -> a)
    -> System m n
mkSystem m f u = System
                    -- < convert from      | actual thing | convert to >
    { sysInertia       =                     m
    , sysCoords        =        gvecR      . f            . grVec
    , sysJacobian      = tr   . vec2l      . jacobianT f  . grVec
    , sysHessian       = tr2  . fmap vec2l . hessianF f   . grVec
    , sysPotential     =                     u            . grVec
    , sysPotentialGrad =        gvecR      . grad u       . grVec
    }

-- | Equations of motion for a system at a given position in phase space
hamilEqns
    :: (KnownNat n, KnownNat m)
    => System m n
    -> Phase n
    -> (R n, R n)       -- dq/dt and dp/dt
hamilEqns s (Phase q p) = (dqdt, dpdt)
  where
    j       = sysJacobian s q
    trj     = tr j
    mHat    = diag (sysInertia s)
    kHat    = trj `mul` mHat `mul` j
    kHatInv = inv kHat
    dqdt    = kHatInv #> p
    dpdt    = gvecR bigUglyThing - sysPotentialGrad s q
      where
        bigUglyThing =
          fmap (\j2 -> -p <.> kHatInv #> trj #> mHat #> j2 #> kHatInv #> p)
               (sysHessian s q)

-- | Step a system's position through phase space using Euler integration
stepEuler
    :: (KnownNat n, KnownNat m)
    => System m n       -- ^ the system
    -> Double           -- ^ dt
    -> Phase n          -- ^ q(t) and p(t)
    -> Phase n          -- ^ q(t + dt) and p(t + dt)
stepEuler s dt ph@(Phase q p) = Phase (q + konst dt * dq) (p + konst dt * dp)
  where
    (dq, dp) = hamilEqns s ph

-- | Iterate Euler's method to repeatedly step a system through phase space
runSystem
    :: (KnownNat n, KnownNat m)
    => System m n       -- ^ the system
    -> Double           -- ^ dt
    -> Phase n          -- ^ initial phase
    -> [Phase n]        -- ^ progression of the system using Euler integration
runSystem s dt = go
  where
    go p0 = p0 : go (stepEuler s dt p0)

-- | A simple particle in 2D cartesian space under gravity
simpleSystem :: System 2 2
simpleSystem = mkSystem (vec2 5 5) id pot
  where
    -- potential energy of a gravity field
    -- U(x,y) = 9.8 * y
    pot :: RealFloat a => V.Vector 2 a -> a
    pot xy = 9.8 * (xy `V.index` 1)

-- | An initial position in configuration space, representing a particle
-- at <0,0> with initial velocity <1,3>
simpleConfig0 :: Config 2
simpleConfig0 = Config
    { confPositions  = vec2 0 0
    , confVelocities = vec2 1 3
    }

simpleMain :: IO ()
simpleMain =
    mapM_ (disp 2 . phasePositions)  -- position with 2 digits of precision
  . take 25                          -- 25 steps
  $ runSystem simpleSystem 0.1 (toPhase simpleSystem simpleConfig0)

-- | A pendulum system, parameterized by its angle clockwise from
-- equilibrium
pendulum :: System 2 1
pendulum = mkSystem (vec2 5 5) coords pot      -- 5kg particle
  where
    -- <x,y> = <-0.5 sin(theta), -0.5 cos(theta)>
    -- pendulum of length 0.25
    coords :: RealFloat a => V.Vector 1 a -> V.Vector 2 a
    coords (V.head->theta) = V.fromTuple (- 0.25 * sin theta, - 0.25 * cos theta)
    -- potential energy of gravity field
    -- U(x,y) = 9.8 * y
    pot :: RealFloat a => V.Vector 1 a -> a
    pot q = 9.8 * (coords q `V.index` 1)

-- | An initial pendulum position in configuration space, representing
-- a pendulum at equilibrium with initial angular velocity 0.1 rad/sec
-- clockwise
pendulumConfig0 :: Config 1
pendulumConfig0 = Config
    { confPositions  = 0
    , confVelocities = 0.1
    }

pendulumMain :: IO ()
pendulumMain =
    mapM_ (disp 3 . phasePositions)  -- position with 2 digits of precision
  . take 25                          -- 25 steps
  $ runSystem pendulum 0.1 (toPhase pendulum pendulumConfig0)

main :: IO ()
main = pendulumMain

max_stars_repo_name: mstksg/inCode

max_issues_repo_name: mstksg/inCode

max_forks_repo_name: mstksg/inCode


text:
{-# LANGUAGE GADTs #-}
module Math.IRT.Model.FourPLM
  ( FourPLM (..)
  ) where

import Numeric.AD (Mode, Scalar, auto)
import Numeric.AD.Mode.Forward.Double
import qualified Numeric.AD.Mode.Tower as T

import Statistics.Distribution

import Math.IRT.Internal.Distribution
import Math.IRT.Internal.LogLikelihood
import Math.IRT.Model.Generic

data FourPLM = FourPLM { discrimination :: !Double
                       , difficulty     :: !Double
                       , pseudoGuessing :: !Double
                       , asymptote      :: !Double
                       } deriving (Show)

instance Distribution FourPLM where
    cumulative = cumulative4PL

instance ContDistr FourPLM where
    density  x = diff (cumulative4PL x)
    quantile _ = error "This shouldn't be needed"

instance DensityDeriv FourPLM where
    densityDeriv x = (!! 2) . T.diffs (cumulative4PL x)

instance GenericModel FourPLM where
    fromRasch         b = FourPLM 1.0 b 0.0 1.0
    fromOnePLM        b = FourPLM 1.7 b 0.0 1.0
    fromTwoPLM      a b = FourPLM   a b 0.0 1.0
    fromThreePLM  a b c = FourPLM   a b   c 1.0
    fromFourPLM         = FourPLM

instance LogLikelihood FourPLM where
    logLikelihood = logLikeFunc cumulative4PL

cumulative4PL :: (Mode a, Floating a, Scalar a ~ Double) =>
       FourPLM -- ^The IRT parameters
    -> a       -- ^Theta
    -> a
cumulative4PL params theta =
    let (a, b, c, d) = makeParamAutos params in
    c + (              (d - c)
         -----------------------------------
         / (1 + exp ((-a) * (theta - b))))
    where makeParamAutos (FourPLM sa sb sc sd) =
              let a = auto sa
                  b = auto sb
                  c = auto sc
                  d = auto sd
              in (a, b, c, d)

max_stars_repo_name: argiopetech/irt

max_issues_repo_name: argiopetech/irt

max_forks_repo_name: argiopetech/irt


text:
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE DeriveGeneric #-}
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeOperators #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-#LANGUAGE OverloadedLists#-}

module Lib
    ( app
    , main'
    ) where

import Network.Wai
import Network.Wai.Handler.Warp
import Servant
import Servant.API
import Data.Aeson
import GHC.Generics
import Data.String.Conversions
import Text.Read (readMaybe)
import Servant.HTML.Lucid
import Lucid
import Data.Time.Calendar
import Data.List (sortBy)
import Data.Complex

import Data.HashMap.Strict
import Data.Text(Text(..))

------------------------ API SPEC ---------------------------
type UserAPI =
           "users" :> QueryParam "sortby" SortBy :> Get '[JSON] [User]
      :<|> "members" :> Get '[JSON] (Headers '[Header "Access-Control-Allow-Origin" String] Members)
      :<|> "albert" :> Get '[JSON] User
      :<|> "isaac" :> Get '[JSON] User

type MathAPI =
          "sum" :> Capture "x" Int :> Capture "y" Int :> Get '[JSON] Int
     :<|> "conjugate" :> Capture "x" (Complex Int)  :> Get '[JSON] String

type MyAPI = UserAPI :<|> MathAPI

data User = User
    { name :: String
    , age :: Int
    , img :: String
    } deriving (Eq, Show, Generic)

newtype Members = Members [User]

--------------------- INSTANCES ---------------------
instance ToJSON User

instance ToJSON Members where
   toJSON (Members us) = Object [("members", toJSON us)]

instance ToHtml User where
    toHtml u = tr_ $ do
                td_ $ toHtml (name u)
                td_ $ toHtml (show $ age u)

    toHtmlRaw = toHtml

instance ToHtml [User] where
    toHtml u = table_ $ do
        tr_ $ do
            th_ $ toHtml ("Name"::String)
            th_ $ toHtml ("Age"::String)
        foldMap toHtml u

    toHtmlRaw = toHtml

data SortBy = Age | Name deriving (Generic, Show,Eq)

instance FromHttpApiData SortBy where
    parseQueryParam "age" = Right Age
    parseQueryParam "name" = Right Name
    parseQueryParam _ = Left "Not valid query param"

instance (Num a,Read a) => FromHttpApiData (Complex a) where
    --parseUrlPiece :: Text -> Either Text a
    parseUrlPiece txt = maybe (Left "Cant parse complex num") (Right) (readMaybe $ cs txt)

---------------------- ENDPOINT HANDLERS -------------------
users :: [User]
users =
  [ User "Isaac Newton"    372 "http://placehold.it/250x250"
  , User "Albert Einstein" 136 "http://placehold.it/250x250"
  , User "Haskell Curry" 200 "http://placehold.it/250x250"
  ]

members :: Monad m => m (Headers '[Header "Access-Control-Allow-Origin" String] Members)
members = return $ addHeader "*" (Members users)

albert = users !! 0
isaac = users !! 1

sortedUsers Age = sortBy (compareWith age)
sortedUsers Name = sortBy (compareWith name)

-- m is actually a EitherT ServantErr IO, in the abscense of SortBy query param users gets returned
sortUsers :: Monad m =>  Maybe SortBy -> m [User]
sortUsers = return . maybe users (flip sortedUsers users)

server1 :: Server UserAPI
server1 =
          sortUsers
     :<|> members
     :<|> return isaac
     :<|> return albert

server2 :: Server MathAPI
server2 = sum
      :<|> showConjugate
      where
       sum = return <.. (+)

       showConjugate :: Monad m => Complex Int -> m String
       showConjugate = return . show . conjugate

server :: Server MyAPI
server = server1 :<|> server2

myAPI :: Proxy MyAPI
myAPI = Proxy

---------------------------------- SERVER --------------------------------
{- 'serve' comes from servant and hands you a WAI Application,
which you can think of as an "abstract" web application,
not yet a webserver. -}

app :: Application
app = serve myAPI server

main' = do
    let hostIP = "0.0.0.0"
    let port = 5000
    putStrLn $ "Server running in " ++ (show hostIP) ++ " port " ++ (show port)
    let settings = setPort port $ setHost hostIP defaultSettings
    runSettings settings app

----------------------------- UTILITIES -----------------------------
(<..) = (.) . (.)
compareWith f x y = compare (f x) (f y)

max_stars_repo_name: DanielCardonaRojas/LearnServant

max_issues_repo_name: DanielCardonaRojas/LearnServant

max_forks_repo_name: DanielCardonaRojas/LearnServant


text:
{-# OPTIONS_GHC -fplugin Numeric.Units.Dimensional.DK.Solver #-}

{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE TypeOperators #-}

module Controls
where

import Numeric.Units.Dimensional.DK.Prelude
import Numeric.LinearAlgebra.Dimensional.DK.Shapes
import qualified Prelude as P

data DimMat (s :: Shape) (e :: *) = DimMat

(<>) :: DimMat s1 e -> DimMat s2 e -> DimMat (ShapeProduct s1 s2) e
_ <> _ = DimMat

add :: DimMat s e -> DimMat s e -> DimMat s e
add _ _ = DimMat

data ContinuousLiSystem (iv :: Dimension) (xs :: [Dimension]) (ys :: [Dimension]) (us :: [Dimension]) e = ContinuousLiSystem
                                                                                {
                                                                                a'' :: DimMat (DivideVectorLists (MapDiv iv xs) xs) e,
                                                                                b'' :: DimMat (DivideVectorLists (MapDiv iv xs) us) e,
                                                                                c'' :: DimMat (DivideVectorLists ys xs) e,
                                                                                d'' :: DimMat (DivideVectorLists ys us) e
                                                                                }

type ContinuousLtiSystem = ContinuousLiSystem DTime

type ExampleSystem = ContinuousLtiSystem
    '[DLength, DVelocity, DPlaneAngle, DAngularVelocity]
    '[DLength, DPlaneAngle]
    '[DForce]
    Double

evaluate :: ContinuousLiSystem iv (x ': xs) (y ': ys) (u ': us) e -> DimMat (VectorShape x xs) e -> DimMat (VectorShape u us) e -> (DimMat (VectorShape (x / iv) (MapDiv iv xs)) e,  DimMat (VectorShape y ys) e)
evaluate sys x u = let
                      a = a'' sys
                      b = b'' sys
                      c = c'' sys
                      d = d'' sys
                      xDot = (a <> x) `add` (b <> u)
                      y = (c <> x) `add` (d <> u)
                    in (xDot, y)

{- Example from http://ctms.engin.umich.edu/CTMS/index.php?example=InvertedPendulum&section=ControlStateSpace -}

-- not really stated on that page, but if you go back a couple of pages in their derivation
-- you can see that the type of u is a 1x1 matrix whose sole element is a force

massOfCart = (0.5 :: Double) *~ (kilo gram)
massOfPendulum = (0.2 :: Double) *~ (kilo gram)
coefficientOfFrictionForCart = (0.1 :: Double) *~ (newton / (meter / second))
lengthToPendulumCenterOfMass = (0.3 :: Double) *~ meter
massMomentOfInertiaOfPendulum = (0.006 :: Double) *~ (kilo gram * meter^pos2)
g = (9.8 :: Double) *~ (meter / second^pos2)

p = massMomentOfInertiaOfPendulum*(massOfCart+massOfPendulum)+(massOfCart*massOfPendulum*lengthToPendulumCenterOfMass*lengthToPendulumCenterOfMass)

a22 = negate (massMomentOfInertiaOfPendulum+massOfPendulum * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass) * coefficientOfFrictionForCart / p
a23 = (massOfPendulum * massOfPendulum * g * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass) / p
a42 = negate (massOfPendulum * lengthToPendulumCenterOfMass * coefficientOfFrictionForCart) / p
a43 = massOfPendulum * g * lengthToPendulumCenterOfMass*(massOfCart + massOfPendulum)/p

b21 = (massMomentOfInertiaOfPendulum + (massOfPendulum * lengthToPendulumCenterOfMass * lengthToPendulumCenterOfMass)) / p
b41 = massOfPendulum * lengthToPendulumCenterOfMass / p

{-
pendulum = ContinuousLiSystem {
           a'' = [matD| _0, _1, _0, _0;
                        _0, a22, a23, _0;
                        _0, _0, _0, _1;
                        _0, a42, a43, _0 |],
           b'' = [matD| _0;
                        b21;
                        _0;
                        b41 |],
           c'' = [matD| _1, _0, _0, _0;
                        _0, _0, _1, _0 |],
           d'' = zeroes
         } :: ExampleSystem
-}

pendulum = ContinuousLiSystem {
           a'' = DimMat,
           b'' = DimMat,
           c'' = DimMat,
           d'' = DimMat
         } :: ExampleSystem

max_stars_repo_name: dmcclean/dimensional-dk-linalg

max_issues_repo_name: dmcclean/dimensional-dk-linalg

max_forks_repo_name: dmcclean/dimensional-dk-linalg


text:
{-# LANGUAGE ScopedTypeVariables, ConstraintKinds, GADTs #-}
{-# LANGUAGE FunctionalDependencies #-}
{-# LANGUAGE UndecidableInstances #-} -- See below
{-# LANGUAGE TupleSections #-}
{-# OPTIONS_GHC -Wall #-}
{-# OPTIONS_GHC -fno-warn-unused-imports #-}
{-# OPTIONS_GHC -fno-warn-unused-binds   #-}

----------------------------------------------------------------------
-- |
-- Module      :  Generic.FFT
-- Copyright   :  (c) 2013 Tabula, Inc.
--
-- Maintainer  :  conal@tabula.com
-- Stability   :  experimental
--
-- FFT via functor composition. Apply to depth-typed perfect binary leaf trees.
----------------------------------------------------------------------
module Generic.FFT
  ( dft
  , idft
  , dct
  , HasFFT(..)
  ) where

import Prelude hiding (sum)

import Control.Applicative (Applicative(..), liftA2)
import Control.Arrow ((***))
import Data.Foldable (Foldable, sum)
import Data.Functor ((<$>))
import Data.Monoid (Monoid(..), Sum(..), (<>))
import Data.Traversable (Traversable(..), mapAccumL)

import Data.Complex (Complex(..), conjugate, realPart)
import Data.Tuple (swap)

import Data.UniformPair

import TypeUnary.Nat

import qualified Data.FTree.BottomUp as B
import qualified Data.FTree.TopDown as T

{--------------------------------------------------------------------
    Misc
--------------------------------------------------------------------}
type Unop a = a -> a

transpose :: (Traversable g, Applicative f) => g (f a) -> f (g a)
transpose = sequenceA

inTranspose ::
     (Traversable f, Traversable k, Applicative g, Applicative h)
  => (g (f a) -> k (h b))
  -> (f (g a) -> h (k b))
inTranspose = transpose --> transpose

infixr 1 -->

-- | Add pre- and post processing
(-->) :: (a' -> a) -> (b -> b') -> ((a -> b) -> (a' -> b'))
(f --> h) g = h . g . f

type R = Double
type C = Complex R

scanL :: (Traversable f, Monoid a) => (a, f a) -> (f a, a)
scanL (a0, as) = swap (mapAccumL h a0 as)
  where
    h a a' = (a <> a', a)

-- TODO: Replace scanL with my efficient parallel version.
-- Prefix (left) sums
sumsL :: (Traversable f, Num a) => (a, f a) -> (f a, a)
sumsL = (fmap getSum *** getSum) . scanL . (Sum *** fmap Sum)

-- Yield a structure counting from 0 to size-1, together with size
counts ::
     forall f a. (TA f, Num a)
  => (f a, a)
counts = sumsL (0, pure 1)

-- Cross product of structures
cross :: (Functor g, Functor f) => g a -> f b -> g (f (a, b))
as `cross` bs = fmap (\a -> fmap (a, ) bs) as

-- All products of numbers from each structure.
products :: (Functor g, Functor f, Num n) => g n -> f n -> g (f n)
products = (fmap . fmap . fmap . fmap) (uncurry (*)) cross

-- as `products` bs = (fmap.fmap) (uncurry (*)) (as `cross` bs)
-- as `products` bs = fmap (\ a -> fmap (\ b -> a*b) bs) as
-- Dot product of structures. Assumes a trie-like f, so that the Applicative
-- instance combines corresponding elements. Perhaps replace Applicative with a
-- more suitable constraint.
dot :: (Applicative f, Foldable f, Num a) => f a -> f a -> a
u `dot` v = sum (liftA2 (*) u v)

i2pi :: C
i2pi = 0 :+ 2 * pi

-- Principle nth root of unity with negated angle
uroot :: Int -> C
uroot n = exp (-i2pi / fromIntegral n)

{--------------------------------------------------------------------
    Unoptimized discrete Fourier transform (DFT)
--------------------------------------------------------------------}
-- Discrete Fourier transform:
-- $X_k = \sum_{n=0}^{N-1} x_n e^{-i 2\pi k \frac{n}{N}}$ for $k = 0,\ldots,N$.
dft :: TA f => Unop (f C)
dft xs = (xs `dot`) <$> rootses

idft ::
     forall f. TA f
  => Unop (f C)
idft = fmap ((/ fromIntegral n) . conjugate) . dft . fmap conjugate
  where
    n :: Int
    indices :: f Int
    (indices, n) = counts

-- Inefficinet; doesn't exploit symmetry, but a starting point for testing
dct :: TA f => Unop (f R)
dct = fmap realPart . dft . fmap (:+ 0)

-- Powers of 'uroot' needed in the DFT:
-- $e^{\frac{-i 2\pi k n}{N}}$ for $k,n = 0,\ldots,N$:
rootses ::
     forall f. TA f
  => f (f C)
rootses = rootCross tot indices indices
  where
    indices :: f Int
    (indices, tot) = counts

rootCross :: (Functor g, Functor f, Integral n) => Int -> g n -> f n -> g (f C)
rootCross tot = (fmap . fmap . fmap . fmap) (uroot tot ^) products

-- rootCross tot is js = (fmap.fmap) (uroot tot ^) (is `products` js)
{--------------------------------------------------------------------
    FFT
--------------------------------------------------------------------}
-- | FFT computation, parametrized by structure
class HasFFT f f' | f -> f' where
  fft :: f C -> f' C

-- Constraint shorthands
type TA f = (Traversable f, Applicative f)

type TAH f f' = (TA f, TA f', HasFFT f f')

-- Binary butterfly.
instance HasFFT Pair Pair where
  fft (a :# b) = a + b :# a - b

-- Decimation in time (DIT)
instance (TAH f f', IsNat n) => HasFFT (B.T f n) (T.T f' n) where
  fft (B.L a) = T.L a
  fft (B.B t) = T.B (fftC t)

-- Decimation in frequency (DIF)? I'm unsure.
instance (TAH f f', IsNat n) => HasFFT (T.T f n) (B.T f' n) where
  fft (T.L a) = B.L a
  fft (T.B t) = B.B (fftC t)

instance HasFFT [] [] where
  fft = dft

--     Variable s `f, f' occur more often than in the instance head
--       in the constraint: TAH f
--     (Use -XUndecidableInstances to permit this)
--     In the instance declaration for `HasFFT (T f n)'
--
-- This warning vanishes when we spell out TAH. Hm.
-- I'd prefer terser definitions like `fft = T.inT' T.l (T.B . fftC)`, but I
-- haven't found a type for `inT'` that GHC likes.
fftsT :: (Applicative f, Traversable g, HasFFT g g') => g (f C) -> f (g' C)
fftsT = fmap fft . transpose

-- FFT of composed functors
fftC :: (TAH f f', TAH g g') => g (f C) -> f' (g' C)
fftC = transpose . fftsT . twiddle . fftsT

{-

--   fftsT     :: g  (f  C) -> f  (g' C)
--   twiddle   :: f  (g' C) -> f  (g' C)
--   fftsT     :: f  (g' C) -> g' (f' C)
--   transpose :: g' (f' C) -> f' (g' C)

-}
-- Multiply by twiddle factors
twiddle ::
     forall f g. (TA f, TA g)
  => Unop (g (f C))
twiddle = (liftA2 . liftA2) (*) rootses'
  where
    rootses' = rootCross (gTot * fTot) gIndices fIndices
      where
        fIndices :: f Int
        (fIndices, fTot) = counts
        gIndices :: g Int
        (gIndices, gTot) = counts

{--------------------------------------------------------------------
    Experimental variation
--------------------------------------------------------------------}
-- This version is a better fit with
-- <https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm#General_factorizations>
-- FFT of composed functors
fftC' :: (TAH f f', TAH g g') => g (f C) -> f' (g' C)
fftC' = fftsT' . transpose . twiddle . fftsT'
  where
    fftsT' :: (TA h, TAH k k') => k (h C) -> k' (h C)
    fftsT' = (inTranspose . fmap) fft

-- Hm. This definition differs from the previous one, since `twiddle` and
-- `transpose` got swapped. I doubt they're equivalent.
-- Types in the fftC' definition (right to left):
--
--   fftsT'    :: g  (f  C) -> g' (f  C)
--   twiddle   :: g' (f  C) -> g' (f  C)
--   transpose :: g' (f  C) -> f  (g' C)
--   fftsT'    :: f  (g' C) -> g' (g' C)
{--------------------------------------------------------------------
    Tests
--------------------------------------------------------------------}
-- {1, 0, 0, 0, ...} <=DFT=> {1, 1, 1, 1, ...}\
-- {1, -1, 1, -1, 1, -1, ...} <=DFT=> {0, 0, ... , N, 0, 0, ...} (where the 'N' occurs at the N/2 position in the output vector.
-- dft . idft ~= id
-- TODO :: show some example outputs and add a series of Quickcheck tests to
-- validate them.
-- TODO: make B / Pair showable

_p1 :: Pair C
_p1 = 1 :# 0

_t1 :: B.T Pair N3 C
_t1 = B.B (B.B (B.B (B.L (((1 :# 0) :# (0 :# 0)) :# ((0 :# 0) :# (0 :# 0))))))

_t2 :: B.T Pair N4 C
_t2 = pure 1

_t3 :: B.T Pair N3 C
_t3 = B.B (pure (1 :# -1))