text:
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE BangPatterns #-}
-- | Fourier amplitude sensitivity test, a quasi-monte carlo method
--
-- does not yet play well with "ErrorProp.MonteCarlo", though a function
-- like 'ErrorProp.MonteCarlo.sample' should be straightforward to
-- implement. However, 'fast' provides the core of such a function
module ErrorProp.QuasiMC where

import qualified Data.Map as M
import Math.Polynomial
import Math.Polynomial.Bernoulli
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector.Generic as VG
import Data.List

import Data.Numbers.Primes
import Data.VectorSpace

import ErrorProp.QuasiMC.Cukier (cukierOmegas)

import Statistics.Sample
import Data.Maybe

import ErrorProp.MonteCarlo
import ErrorProp.Common

f1 alpha2 = constPoly 1 ^+^ ((-1)^(1 + alpha2 `mod` 2) * (2*pi)^(2*alpha2)) *^ bx
  where bx = bernoulliPoly !! (2*alpha2)

-- the explicit equations for F2,F4,F6 (equation 4.14,15,16) are correct
prop_f11 x = f11 x ~= f1 1 `evalPoly` x
  where f11 x = 1 + 2 * pi^2 * (x^2 - x + 1/6)

prop_f12 x = f12 x ~= f1 2 `evalPoly` x
  where f12 x = 1 + pi^4/45 * (1 - 30* x^2 *  (1-x)^2)

prop_f13 x = f13 x ~= f1 3 `evalPoly` x
  where f13 x = 1 + (2*pi^6/945) * (1 - 21*x^2 + 105*x^4 -126*x^5 + 42*x^6)

x ~= y = (2 * abs(x-y) / (1 + abs x + abs y)) < (1e-8 :: Double)
infix 4 ~=

-- eqn 4.10. Is the function integrated to generate the "best"
-- lattice points
f2 alpha2 xs = V.product $ V.map (f1 alpha2 `evalPoly`) xs

-- | Korbov's construction for good lattice points. Is an LCG
--
-- 4.32
korbov l n s = V.iterateN s (\p -> (p * l) `mod` n) 1

-- | try all values of L up to n/2
korbovLExhaustive alpha2 n s =
    head $
    sort [ (abs (r1-1) :: Double, l)
            | l <- [1 .. n `div` 2 ],
              let r1 = rank1 (korbov l n s) n (V.replicate s 0) (f2 alpha2) ]

-- | assume L has k prime factors. Provides candidate L values,
-- and the error in the integral used to decide that L is the best.
-- The error is absolute and relative, since the target function
-- integrates to 1.
korbovLCompositeK
  :: Int -- ^ alpha/2, roughness of the function being integrated
  -> Int -- ^ k
  -> Int -- ^ n, number of points in the quadrature
  -> Int -- ^ s, dimension
  -> [(Double, Int)] -- ^ (error, best l)
korbovLCompositeK alpha2 k n s =
    sort [ (abs (rankL l - 1) :: Double, l)
            | l <- go k ps 1 ]
  where
    rankL l = rank1 (korbov l n s) n (V.replicate s 0) (f2 alpha2)
    go 0 _ n = [n]
    go k ps1 n = do
      q : ps1' @ (_:_) <- tails ps1
      go (k-1) ps1' (n*q)

    ps = takeWhile ( < round ( (fromIntegral n / 2) ** (1/fromIntegral k)))
                primes

{- | rank1 shifted lattice rule

Given

> f :: [0,1)^s -> R

calculate

> Q(Z,N,C)f = (1/N) sum_j=0^N f { j z/N + c }

is an approximation to the integral of @f@ (@If@) over that region.

equation 4.42 Sloan&Joe 1994 pg 89

If C is multivariant uniform [0,1)^s, we can do a CI for the actual
integral

P( |Q - I| < sigma v) >= 1 - (1/v^2)

with sigma^2 estimated by sum_k=1^q (Q - Qbar)^2 / ( q (q-1))

is unbiased. Quantiles are biased (Lemieux 2008). Bootstrap/jacknife resample to estimate bias.

trapezoidal / rectangle rule is good for the periodic integral here,
but possibly try "GA Evans, JR Webster 1999"

-}
rank1 :: (Integral i, RealFrac a, VG.Vector v a, VG.Vector vi i,
            VG.Vector vi a)
  => vi i -- ^ z
  -> Int -- ^ N
  -> v a -- ^ c
  -> (v a -> a) -- ^ f
  -> a
rank1 z n c f = go 0 0
  where
    zn = VG.convert $ VG.map (\zi -> fromIntegral zi / fromIntegral n) z

    go !j !accum
      | j == VG.length zn = accum
      | otherwise = go (j+1) $ f $ VG.zipWith (+) c
              $ VG.map (\zni -> fracPart $ fromIntegral j*zni) zn

{- |

>>> map fracPart [12.3, -2.5]
[0.3, 0.5]

-}
fracPart x = x - fromInteger (floor x)

-- | fourier amplitude sensitivity test
fast ::
  (RealFrac b, VG.Vector V.Vector b) =>
    V.Vector Int -- ^ N incommensurate frequencies (see 'korbov', 'cukierOmegas', or others)
  -> (V.Vector b -> Double) -- ^ a function of N variables which will be sampled on @[0, 1)^N@
  -> V.Vector Double -- ^ fraction of the total variance attributed to the i'th variable
fast omegas f =
    let omegaMax = VG.maximum omegas
        q = omegaMax

        xj k = VG.map (\w -> fracPart $ fromIntegral (w*k) / fromIntegral (2*q+1)) omegas
        sks = V.enumFromTo (-q) q
        fks = VG.map (f . xj) sks

        withK k = VG.map (\w -> mean
                              $ VG.zipWith (\s f -> f * k (w*s)) sks fks)
                omegas

        nq = fromIntegral (2*q+1)

        a i = withK (cos . (*(2*i*pi/nq)) . fromIntegral)
        b i = withK (sin . (*(2*i*pi/nq)) . fromIntegral)

        vjs = foldl1 (V.zipWith (+))
              $ map (V.map (^2))
              [  if odd j then b fij else a fij
                  | j <- [1 .. 2],
                    let fij = fromIntegral j ]

        -- normalizeFac = 2 / variance fks -- XXX the
        -- correct normalization should be some analytical
        -- expression... like above, not the value below
        normalizeFac = recip $ V.sum vjs

    in V.map (*normalizeFac) vjs

fastKi s np k =
  let n = primes !! np
      (_, l) = korbovLExhaustive 1 n s
      omegas = korbov l n s

  in fast omegas (\v -> v V.! k)

fastKis np s = [ (k,fastKi s np k)  | k <- [0 .. s-1]]
{- properties.

> ki k = fast omegas (\x -> x V.! k)

> ki k V.! j == 1 if k==j else 0

-}

-- | `pickNP err s` uses 'fastKis' to pick the parameter np, such that the
-- value that should be 1 is within err of 1, for a problem of dimension s
pickNP err s = head [ np
    | np <- [ 2 .. ],
     all ok (fastKis np s) ]
  where ok (i, v) = abs (1 - v V.! i) < err

sampleQMC :: Double -- ^ error tolerance (see 'pickNP')
  -> MV Double
  -> [(String, Double)] -- variable name, fraction of the total variance
sampleQMC err mv =
  let n = wheelSieve 5 !! pickNP err s
      (_, l) = korbovLExhaustive 1 n s
      omegas = korbov l n s
      s = M.size (mv_vars mv)

      f v = mv_sample mv (updateMV v (mv_vars mv))
  in M.keys (mv_vars mv) `zip` VG.toList (fast omegas f)

updateMV :: (Floating b, V.Unbox b) => V.Vector b -> M.Map String (Maybe b) -> M.Map String b
updateMV v m = snd $ M.mapAccum (\ n e -> (n+1, fromMaybe 1 e * toUnitVar (v V.! n))) 0 m
  where
        toUnitVar x = sqrt 12 * (x - 0.5)

max_stars_repo_name: aavogt/ErrorProp

max_issues_repo_name: aavogt/ErrorProp

max_forks_repo_name: aavogt/ErrorProp


text:
import Test.Hspec
import Numeric.LinearAlgebra

import Statistics.GLM

main :: IO ()
main = hspec $ do
  describe "logReg" $ do
    it "can perform logistic regression" $ do
      fmap (truncate . (*10)) (toList (last $ take 5 $ glm logisticReg testA testX testB)) `shouldBe` fmap (truncate . (*10)) (toList $ vector [log 8, log 0.25])

testA :: Matrix Double
testA = (6><2) [1::Double, 1, 1, 1, 1, 1, 1, 2, 1, 2, 1, 2]

testB :: Vector Double
testB = vector [1, 1, 0, 1, 0, 0]

testX :: Vector Double
testX = vector [0, 0]

ext: hs


text:
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# LANGUAGE UndecidableInstances #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE TypeFamilies #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}

module HLearn.Models.Distributions.Univariate.Binomial
    where

import HLearn.Algebra
import HLearn.Models.Distributions.Common
import HLearn.Models.Distributions.Univariate.Internal.Moments
import HLearn.Models.Distributions.Visualization.Gnuplot

import qualified Statistics.Distribution as S
import Statistics.Distribution.Binomial

-------------------------------------------------------------------------------
-- data types

newtype Binomial sample prob = Binomial {  bmoments :: (Moments3 sample) }
    deriving (Read,Show,Eq,Ord,Monoid,Group)

-------------------------------------------------------------------------------
-- Training

instance (Num sample) => HomTrainer (Binomial sample prob) where
    type Datapoint (Binomial sample prob) = sample
    train1dp dp = Binomial $ train1dp dp

-------------------------------------------------------------------------------
-- distribution

instance (Num sample) => Probabilistic (Binomial sample prob) where
    type Probability (Binomial sample prob) = prob

instance (Floating prob) => PDF (Binomial Int Double) where
    pdf (Binomial dist) dp = S.probability (binomial n p) dp
        where
            n = bin_n $ Binomial dist
            p = bin_p $ Binomial dist

bin_n :: Binomial Int Double -> Int
bin_n (Binomial dist) = round $ ((fromIntegral $ m1 dist :: Double) / (fromIntegral $ m0 dist)) / (bin_p $ Binomial dist)

bin_p :: Binomial Int Double -> Double
bin_p (Binomial dist) = ((fromIntegral $ m1 dist) / (fromIntegral $ m0 dist)) + 1 - (fromIntegral $ m2 dist)/(fromIntegral $ m1 dist)

instance
    ( PDF (Binomial sample prob)
--     , PlottableDataPoint sample
    , Show prob
    , Show sample
    , Ord sample
    , Ord prob
    , Num prob
    , Integral sample
    ) => PlottableDistribution (Binomial sample prob)
-- instance PlottableDistribution (Poisson Int Double)
        where

    plotType _ = Points

    samplePoints dist = [min..max]
        where
            min = 0
            max = maximum [20,floor $ 3*mu]
            mu = 5

-- instance (Fractional prob) => Mean (Binomial prob) where
--     mean (Binomial dist) =
--
-- instance (Fractional prob) => Variance (Binomiral prob) where
--     variance dist = mean dist


text:
import Control.Monad (liftM)
import Control.Parallel
import Data.Binary.Get (runGet)
import Graphics.Rendering.Cairo hiding (Matrix)
import Numeric.LinearAlgebra
import System.Random (mkStdGen, randoms)
import qualified Data.ByteString.Lazy as BS
import Gradient
import BinaryExtras
import RandomExtras
import Types
import Unsafe.Coerce
import GHC.Float
import Data.Binary.Get
import Data.Int
import Data.Word

opts = morrowind

main = do
  contents <- BS.readFile $ inputFile opts
  let arr0 = runGet (getManyF (rawFunction opts) $ imageWidth opts*imageHeight opts) contents
  let landscape = (imageHeight opts><imageWidth opts) $ arr0
--  landscape <- liftM ((imageHeight opts><imageWidth opts) . runGet (getAll $ imageWidth opts*imageHeight opts)) $ BS.readFile $ inputFile opts
  let jitter = randoms $ mkStdGen 53
  let [miny, maxy, minx, maxx] = map fromIntegral [1, imageHeight opts-2, 1, imageWidth opts-2]
  let (numpointsy, numpointsx) = (floor $ (maxy-miny)/(gridWidth opts), floor $ (maxx-minx)/(gridWidth opts))
  let regulargrid = [(y,x) | y <- [0..numpointsy-1], x<-[0..numpointsx-1]]
  let jitteredgrid = zipWith (\(y,x) (jy,jx) ->
        ( miny + (gridWidth opts)* (fromIntegral y + jy)
        , minx + (gridWidth opts) * (fromIntegral x + jx)
        )) regulargrid jitter
  let ascents = map (gradientStepper Ascent (ascentDelta opts) landscape) jitteredgrid
  s <- createImageSurface FormatARGB32
    (floor $ multiplier opts * (fromIntegral $ imageWidth opts))
    (floor $ multiplier opts * (fromIntegral $ imageHeight opts))
  renderWith s $ do
    setSourceRGBA 1 1 1 0.0
    paint
  mapM_ (drawLine (lineWidth opts) (cellWidth opts) (multiplier opts) s) ascents
  surfaceWriteToPNG s $ outputFile opts

drawLine line cell mult s lista = do
  renderWith s $ do
    setLineWidth line
    setLineCap LineCapRound
    go lista
  where
  go ((x1,y1,z1):loput@((x2,y2,z2):_)) = do
    let c = 1.0 - ((atan2 (abs(z1-z2)) cell)/(pi/2)) in setSourceRGB c c c
    moveTo (x1*mult) (y1*mult)
    lineTo (x2*mult) (y2*mult)
    stroke
    go loput
  go _ = return ()

-- Example usage

morrowind = Options
  { inputFile = "tesannwyn.raw"
  , imageWidth = 2688
  , imageHeight = 2816
  , gridWidth = 2.5
  , multiplier = 1.0
  , ascentDelta = 1.0
  , outputFile = "vv.png"
  , lineWidth = 0.35
  , cellWidth = 8.0
  , rawFunction = liftM ((fromIntegral :: Int16 -> Double) . (fromIntegral :: Word16 -> Int16)) getWord16le
  }

oblivion = Options
  { inputFile = "oblivion.raw"
  , imageWidth = 4288
  , imageHeight = 4128
  , gridWidth = 5
  , multiplier = 1.0
  , ascentDelta = 2.0
  , outputFile = "obl.png"
  , lineWidth = 0.5
  , cellWidth = 8.0
  , rawFunction = liftM ((fromIntegral :: Int16 -> Double) . (fromIntegral :: Word16 -> Int16)) getWord16le
  }

seworld = Options
  { inputFile = "seworld.raw"
  , imageWidth = 2048
  , imageHeight = 2048
  , gridWidth = 2.5
  , multiplier = 1.0
  , ascentDelta = 1.0
  , outputFile = "se.png"
  , lineWidth = 0.3
  , cellWidth = 8.0
  , rawFunction = liftM ((fromIntegral :: Int16 -> Double) . (fromIntegral :: Word16 -> Int16)) getWord16le
  }

mokki = Options
  { inputFile = "2mdem.flt"
  , imageWidth = 1857
  , imageHeight = 2136
  , gridWidth = 12.5
  , multiplier = 1.5
  , ascentDelta = 2.0
  , outputFile = "mok.png"
  , lineWidth = 0.5
  , cellWidth = 0.125
  , rawFunction = liftM (float2Double . unsafeCoerce) getWord32le
  }

turku = Options
  { inputFile = "pohja"
  , imageWidth = 513
  , imageHeight = 513
  , gridWidth = 1.25
  , multiplier = 2
  , ascentDelta = 1.0
  , outputFile = "turku.png"
  , lineWidth = 0.5
  , cellWidth = 8
  , rawFunction = liftM ((fromIntegral :: Int16 -> Double) . (fromIntegral :: Word16 -> Int16)) getWord16le
  }

max_stars_repo_name: angs/Hachure-map

max_issues_repo_name: angs/Hachure-map

max_forks_repo_name: angs/Hachure-map


text:
{-# LANGUAGE GADTs #-}
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE TypeApplications #-}
{-# LANGUAGE ViewPatterns #-}
module Main where

import Config.Schema (SectionsSpec, ValueSpec)
import Control.Concurrent
import Control.Concurrent.Chan
import Control.Concurrent.STM
import Control.Concurrent.STM.TVar
import Control.DeepSeq
import Control.Exception
import Control.Monad
import Control.Monad.IO.Class
import Data.Bits
import Data.ByteString (ByteString)
import Data.Char
import Data.Coerce
import Data.Colour.CIE
import Data.Colour.CIE.Illuminant
import Data.Colour.SRGB
import Data.Containers.ListUtils
import Data.Default
import Data.Foldable
import Data.Functor.Alt ((<!>))
import Data.String
import Data.Tree
import Data.Int
import Data.IORef
import Data.List
import Data.Map (Map)
import Data.Maybe
import Data.Monoid
import Data.Multiset (Multiset)
import Data.Ord
import Data.Set (Set)
import Data.Set.Ordered (OSet)
import Data.Text (Text)
import Data.Time
import Data.Time.Calendar.OrdinalDate
import Data.Traversable
import Data.Vector (Vector)
import Data.Void
import Database.PostgreSQL.Simple (Connection)
import GHC.Read
import Graphics.Rendering.Cairo (Render)
import Graphics.UI.Gtk
import Options.Applicative
import System.Clock (Clock(..), TimeSpec(..), getTime)
import System.Directory
import System.Environment
import System.Environment.XDG.BaseDir
import System.FilePath
import System.Glib.Properties
import System.Glib.UTFString
import System.IO
import System.Posix.Process (getProcessID)
import System.Random
import Text.Printf (printf)
import Text.Read
import Text.Read.Lex
import Vertumnion.Shared
import qualified Config as C
import qualified Config.Schema as C
import qualified Data.ByteString as BS
import qualified Data.Map as M
import qualified Data.Map.Strict as M.S
import qualified Data.Multiset as MS
import qualified Data.Set as S
import qualified Data.Set.Ordered as O
import qualified Data.Text as T
import qualified Data.Text.Encoding as T
import qualified Data.Text.IO as T
import qualified Data.Text.Read as T
import qualified Data.Vector as V
import qualified Database.PostgreSQL.Simple as DB
import qualified Database.PostgreSQL.Simple.FromField as DB
import qualified Database.PostgreSQL.Simple.ToField as DB
import qualified Database.PostgreSQL.Simple.FromRow as DB
import qualified Database.PostgreSQL.Simple.ToRow as DB
import qualified Database.PostgreSQL.Simple.Types as DB
import qualified Graphics.Rendering.Cairo as Cairo
import qualified Graphics.Rendering.Cairo.Matrix as CM
import qualified Numeric.LinearAlgebra as Matrix
import qualified Text.PrettyPrint.ANSI.Leijen as Doc

main :: IO ()
main = do
	ctx <- initializeContext
	traverse_ (\thread -> forkIO (thread ctx)) [loggingThread, stdinThread, parserThread, heartbeatThread, moduleThread, preferencesThread]
	guiThread ctx -- NB not forked
	-- TODO: maybe stay open and start ignoring stdin, so that pipes don't fill up and block producers?

data Context = Context
	{ ctxConfig :: Config
	, ctxProfile :: Profile
	, ctxArguments :: Arguments
	, ctxErrors :: Chan String
	, ctxInput :: Chan (TimeSpec, Text)
	-- TODO: since this is only written from parserThread (right?) it can probably be an IORef
	, ctxUITimerLabelStatus :: MVar TimerLabelStatus
	, ctxUITimerLabel :: Label
	, ctxUIGraph :: DrawingArea
	, ctxModuleInput :: Chan Command
	, ctxModuleOutputs :: TVar [(String, TVar [ModulePoint])]
	, ctxTimeMagnitude :: IORef TimeMagnitude
	, ctxEventStore :: TreeStore Event
	, ctxStateRange :: TVar StateRange
	-- the Word is an epoch, used as an equality-check fast path
	, ctxRandomRuns :: TVar (Word, Runs)
	, ctxPreferences :: TChan Preferences
	}

initializeContext :: IO Context
initializeContext = do
	argumentParser <- mkArgumentParser
	arguments <- execParser argumentParser
	profile <- loadProfile (argProfile arguments) (argTargetOverride arguments)
	config <- loadConfig
	errors <- newChan
	input <- newChan
	initGUI
	uiTimerLabelStatus <- newMVar (Stopped 0)
	uiTimerLabel <- labelNew (string <$> Nothing)
	uiGraph <- drawingAreaNew
	moduleInput <- newChan
	moduleOutputs <- newTVarIO []
	timeMagnitude <- newIORef (pDurationHint profile)
	eventStore <- treeStoreNew []
	stateRange <- newTVarIO StateRange
		{ srStateCache = Uninitialized
		, srMin = humanOrder (pTarget profile)
		, srMax = humanOrder (pTarget profile)
		}
	randomRuns <- newTVarIO (0, M.empty)
	preferences <- newTChanIO
	pure Context
		{ ctxConfig = config
		, ctxProfile = profile
		, ctxArguments = arguments
		, ctxErrors = errors
		, ctxInput = input
		, ctxUITimerLabelStatus = uiTimerLabelStatus
		, ctxUITimerLabel = uiTimerLabel
		, ctxUIGraph = uiGraph
		, ctxModuleInput = moduleInput
		, ctxModuleOutputs = moduleOutputs
		, ctxTimeMagnitude = timeMagnitude
		, ctxEventStore = eventStore
		, ctxStateRange = stateRange
		, ctxRandomRuns = randomRuns
		, ctxPreferences = preferences
		}

data Config = Config
	{ db :: ByteString
	} deriving (Eq, Ord, Read, Show)

instance Default Config where
	def = Config
		{ db = "dbname=vertumnion"
		}

configSpec :: ValueSpec Config
configSpec = C.sectionsSpec
	"configuration"
	(config <$> C.optSection "db" help)
	where
	config Nothing = def
	config (Just t) = Config
		{ db = T.encodeUtf8 t -- I know, I know. I don't care, it will never come up
		}
	help = "A postgres connection string for storing splits (default: dbname=vertumnion)"

type ID = Int32

soSort :: SortOrder -> [Text] -> [Text]
soSort = \case
	Ascending -> sortOn humanOrder
	Descending -> sortOn (Down . humanOrder)
	AscendingOn ts -> sortOn (\t -> (O.findIndex t ts, humanOrder t))

soCompare :: SortOrder -> HumanOrder -> HumanOrder -> Ordering
soCompare Ascending t t' = compare t t'
soCompare Descending t t' = compare t' t
soCompare (AscendingOn ts) t t' = case (O.findIndex (toText t) ts, O.findIndex (toText t') ts) of
	(Just i, Just i') -> compare i i'
	_ -> compare t t'

soLt :: SortOrder -> HumanOrder -> HumanOrder -> Bool
soLt so t t' = soCompare so t t' == LT

soMin :: SortOrder -> HumanOrder -> HumanOrder -> HumanOrder
soMin so t t' = if soLt so t t' then t else t'

soMax :: SortOrder -> HumanOrder -> HumanOrder -> HumanOrder
soMax so t t' = if soLt so t t' then t' else t

ctxAllStates :: Context -> Maybe (OSet Text)
ctxAllStates = soAllStates . pSortOrder . ctxProfile

preferencesSpec :: ValueSpec Preferences
preferencesSpec = C.sectionsSpec "preferences" $ pure Preferences
	<*> C.reqSection' "window" rectangleSpec "request that the window manager initially place the window with this geometry"
	<*> C.reqSection' "splits" splitsSpec "set the initial size of the splits pane"
	<*> C.reqSection' "graph" graphSpec "set the initial size of the graph pane"
	where
	rectangleSpec = C.sectionsSpec "window" $ pure Rectangle
		<*> C.reqSection "x" ""
		<*> C.reqSection "y" ""
		<*> C.reqSection "width" ""
		<*> C.reqSection "height" ""
	splitsSpec = C.sectionsSpec "splits" $ C.reqSection "height" ""
	graphSpec = C.sectionsSpec "graph" $ C.reqSection "width" ""

logDir :: IO FilePath
logDir = getUserStateDir "vertumnion"

preferencesDir :: IO FilePath
preferencesDir = configPath "preferences"

loadConfig :: IO Config
loadConfig = do
	fp <- configName
	mh <- try (openFile fp ReadMode)
	case mh of
		Left (e :: IOError) -> def
		Right h -> loadSchema fp h configSpec

loadPreferences :: String -> IO (Maybe Preferences)
loadPreferences profile = do
	fp <- (</> profile) <$> preferencesDir
	mh <- try (openFile fp ReadMode)
	case mh of
		Left (e :: IOError) -> def
		Right h -> Just <$> loadSchema fp h preferencesSpec

data Arguments = Arguments
	{ argProfile :: FilePath
	, argTargetOverride :: Maybe Text
	} deriving (Eq, Ord, Read, Show)

mkArgumentParser :: IO (ParserInfo Arguments)
mkArgumentParser = do
	desc <- usage
	pure $ info (opts <**> helper)
		(  fullDesc
		<> progDesc "Run a splits timer. Events are accepted on stdin."
		<> (footerDoc . Just . Doc.string) desc
		)
	where
	opts = pure Arguments
		<*> strArgument
			(  metavar "PROFILE"
			<> help "A description of the game being played"
			)
		<*> ( optional $ T.pack <$> strOption
				(  long "target"
				<> short 't'
				<> metavar "EVENT"
				<> help "Override the target given in the profile just for this run"
				)
			)

usage :: IO String
usage = do
	prog <- getProgName
	profDir <- profileDir
	prefDir <- preferencesDir
	fp <- configName
	pure . unlines . tail $ [undefined
		, "All files processed by " ++ prog ++ " should be in config-value format and use the"
		, "UTF-8 encoding. See here for details:"
		, "http://hackage.haskell.org/package/config-value/docs/Config.html"
		, ""
		, "PROFILE refers to a file in " ++ profDir ++ "."
		, "The profile should conform to the following format."
		, show (C.generateDocs profileSpec)
		, ""
		, "Preferences will be loaded if they are available. They may be placed in a file"
		, "with the same profile name in " ++ prefDir ++ "."
		, "This file will automatically be created or overwritten during program execution."
		, "The preferences should conform to the following format."
		, show (C.generateDocs preferencesSpec)
		, ""
		, "A configuration file is also read from " ++ fp ++ "."
		, "It is similarly a UTF-8 encoded config-value file, with the following format."
		, show (C.generateDocs configSpec)
		]

data Preferences = Preferences
	{ windowRectangle :: Rectangle
	, splitHeight :: Int
	, graphWidth :: Int
	} deriving (Eq, Show)

scheduleSavePreferences ::
	(WindowClass window, PanedClass graphPane, PanedClass dataPane) =>
	Context -> window -> graphPane -> dataPane -> IO ()
scheduleSavePreferences ctx window graphPane dataPane = liftIO $ do
	(x, y) <- windowGetPosition window
	(w, h) <- windowGetSize window
	sh <- panedGetPosition dataPane
	gw <- panedGetPosition graphPane
	atomically $ writeTChan (ctxPreferences ctx) Preferences
		{ windowRectangle = Rectangle x y w h
		, splitHeight = sh
		, graphWidth = gw
		}

loggingThread :: Context -> IO loop
loggingThread (ctxErrors -> chan) = do
	-- by reading from the channel first, we guarantee not to write a log unless there's something to put in it
	err <- readChan chan
	now <- getCurrentTime
	me <- getProcessID
	dir <- logDir
	let fp = dir </> show me <.> "log"
	mh <- try (createDirectoryIfMissing True dir >> openFile fp AppendMode)
	h <- case mh of
		Right h -> h <$ hSetBuffering h LineBuffering
		Left (_ :: SomeException) -> do
			hPutStrLn stderr $ "WARNING: Could not open log file " ++ fp ++ " for appending. Logging to stderr"
			pure stderr
	hPutStrLn h (show now ++ ": " ++ err)
	forever $ do
		err <- readChan chan
		now <- getCurrentTime
		hPutStrLn h (show now ++ ": " ++ err)

data TimerLabelStatus
	= RunningSince TimeSpec
	| Stopped DiffTimeSpec
	deriving (Eq, Ord, Read, Show)

{-# INLINE string #-}
string :: String -> String
string = id

arbitraryTimeSpec :: TimeSpec
arbitraryTimeSpec = TimeSpec 0 0

guiThread :: Context -> IO ()
guiThread ctx = do
	window <- windowNew
	on window objectDestroy mainQuit
	vbox <- vBoxNew False 0
	set window [windowTitle := string "vertumnion", containerChild := vbox]

	gameLabel <- labelNew . Just . pGame . ctxProfile $ ctx
	-- TODO: Maybe it would be nice for the target to be editable
	targetLabel <- labelNew . Just . ("Target: " <>) . pTarget . ctxProfile $ ctx

	updateTimerLabel ctx

	eventLog <- treeViewNewWithModel (ctxEventStore ctx)
	set eventLog [treeViewHeadersVisible := False]
	let newColumn :: GlibString string => (Event -> IO string) -> IO Int
	    newColumn = addColumnPlainText eventLog (ctxEventStore ctx)
	newColumn (pure . eState)
	newColumn $ \e -> do
		mag <- readIORef (ctxTimeMagnitude ctx)
		pure . snd . showDiffTimeSpec (eMicrosecond e) $ mag
	for_ (pFPS (ctxProfile ctx)) $ \fps -> newColumn $ \e -> do
		mag <- readIORef (ctxTimeMagnitude ctx)
		case eFrame e of
			Nothing -> pure ""
			Just n -> do
				let (mag', content) = showDiffTimeSpec (fromInteger n / realToFrac fps) mag
				when (mag /= mag') $ do
					writeIORef (ctxTimeMagnitude ctx) mag'
					redrawAllIntervals ctx
				pure content

	eventLogScroll <- scrolledWindowNew Nothing Nothing
	set eventLogScroll [containerChild := eventLog, scrolledWindowHscrollbarPolicy := PolicyNever]
	eventLogAdjustment <- scrolledWindowGetVAdjustment eventLogScroll
	-- TODO: disable jumping when position has been manually changed
	onAdjChanged eventLogAdjustment $
		adjustmentGetUpper eventLogAdjustment >>= adjustmentSetValue eventLogAdjustment

	widgetSetSizeRequest (ctxUIGraph ctx) 200 200
	on (ctxUIGraph ctx) draw (renderGraph ctx)

	legend <- mkLegend
	graphPane <- hPanedNew
	panedPack1 graphPane (ctxUIGraph ctx) True True
	panedPack2 graphPane legend False True

	dataPane <- vPanedNew
	panedPack1 dataPane eventLogScroll True True
	panedPack2 dataPane graphPane False True

	boxPackStart vbox gameLabel PackNatural 0
	boxPackStart vbox targetLabel PackNatural 0
	boxPackStart vbox (ctxUITimerLabel ctx) PackNatural 0
	boxPackStart vbox dataPane PackGrow 0

	mprefs <- loadPreferences (argProfile (ctxArguments ctx))
	for_ mprefs $ \prefs -> do
		let Rectangle x y w h = windowRectangle prefs
		windowSetDefaultSize window w h
		windowMove window x y
		panedSetPosition dataPane (splitHeight prefs)
		panedSetPosition graphPane (graphWidth prefs)

	let ssp = scheduleSavePreferences ctx window graphPane dataPane
	on window configureEvent (False <$ liftIO ssp)
	on graphPane (notifyProperty panedPosition') ssp
	on dataPane  (notifyProperty panedPosition') ssp

	widgetShowAll window
	mainGUI

addColumnPlainText ::
	( TreeViewClass treeView
	, TreeModelClass (model row)
	, TypedTreeModelClass model
	, GlibString string
	) =>
	treeView -> model row -> (row -> IO string) -> IO Int
addColumnPlainText treeView model showRow = do
	col <- treeViewColumnNew
	renderer <- cellRendererTextNew
	treeViewColumnPackEnd col renderer False
	cellLayoutSetAttributes col renderer model $ \row -> [cellText :=> showRow row]
	treeViewAppendColumn treeView col

mkLegend :: IO TreeView
mkLegend = do
	store <- listStoreNew $ []
		++ map Left (zip [0..] allModules)
		++ map Right [("time", [1, 1], 1/2), ("state", [1,3,1,1], -1/2), ("number", [1,5], -5/2)]
	treeView <- treeViewNewWithModel store
	col <- treeViewColumnNew
	graphicRenderer <- cellRendererPixbufNew
	labelRenderer <- cellRendererTextNew
	treeViewColumnPackStart col graphicRenderer False
	treeViewColumnPackStart col labelRenderer True
	cellLayoutSetAttributes col graphicRenderer store $ \row -> [cellPixbuf :=> legendPixbuf row]
	cellLayoutSetAttributes col labelRenderer store $ \row -> [cellText := legendText row]
	treeViewAppendColumn treeView col
	set treeView [treeViewHeadersVisible := False]
	pure treeView

legendText :: Either (Double, Module) (String, [Double], Double) -> String
legendText (Left (_, m)) = mLabel m
legendText (Right (s, _, _)) = s

legendPixbuf :: Either (Double, Module) (String, [Double], Double) -> IO Pixbuf
legendPixbuf row = do
	surface <- Cairo.createImageSurface Cairo.FormatARGB32 30 30
	Cairo.renderWith surface $ case row of
		Left (i, _) -> do
			let dAngle = min (pi / numModules) (pi / 12)
			    angle = 2 * pi * i / numModules
			    RGB r g b = toSRGB (cieLAB d65 50 (100 * cos angle) (100 * sin angle))
			Cairo.setSourceRGB r g b
			Cairo.moveTo 15 15
			Cairo.arc 15 15 14 (angle-dAngle) (angle+dAngle)
			Cairo.fill
		Right (_, intervals, start) -> do
			Cairo.scale 5 5
			Cairo.setLineWidth 0.4
			Cairo.setDash intervals start
			Cairo.moveTo 0 3
			Cairo.lineTo 6 3
			Cairo.stroke
	pixbufNewFromSurface surface 0 0 30 30
	where numModules = fromIntegral (length allModules)

renderGraph :: Context -> Render ()
renderGraph ctx = do
	(_, _, w, h) <- Cairo.clipExtents
	let wScale = w/1.75
	    hScale = h/1.3
	if wScale < hScale
		then Cairo.scale wScale wScale >> Cairo.translate 0.45 (h/(2*wScale)-0.55)
		else Cairo.scale hScale hScale >> Cairo.translate (w/(2*hScale)-0.425) 0.1
	states_ <- liftIO (ctxStates ctx)
	ptss_ <- liftIO $ readTVarIO (ctxModuleOutputs ctx)
	ptss <- liftIO $ for ptss_ $ \(label, ptVar) -> (,) label <$> readTVarIO ptVar

	let
		states :: Map Text Double
		states = M.fromList (zip states_ [0..])

		maxStateIx :: Double
		maxStateIx = fromIntegral (M.size states - 1)

		yStates = [s | (_, pts) <- ptss, ModulePoint { y = State s } <- pts]
		yTimes = [t | (_, pts) <- ptss, ModulePoint { y = Time t } <- pts]
		yLogScales = [n | (_, pts) <- ptss, ModulePoint { y = LogScale n } <- pts]

		dtbX = chooseDiffTimeBounds [x pt | (_, pts) <- ptss, pt <- pts]
		dtbY = chooseDiffTimeBounds yTimes
		lsb = chooseLogScaleBounds yLogScales

		isMajor = case pMajorStates (ctxProfile ctx) of
			Nothing -> const True
			Just ss -> (`S.member` ss)

		xPos :: DiffTimeSpec -> Double
		xPos = dtbPos dtbX

		yPos :: Dependent -> Double
		yPos (State s) = 1 - (states M.! s) / maxStateIx
		yPos (Time t) = 1 - dtbPos dtbY t
		yPos (LogScale n) = 1 - logScalePos lsb n

		numModules :: Double
		numModules = fromIntegral (length ptss)

		hasPoints = not . null . concatMap snd $ ptss
		hasState = not . null $ yStates
		hasDiffTime = not . null $ yTimes
		hasLogScale = not . null $ yLogScales

	Cairo.setLineWidth 0.001
	when hasPoints $ do
		let ds = dtbGridLines dtbX

		-- draw grid lines
		for_ ds $ \d -> do
			let coord = xPos d
			Cairo.moveTo coord 0
			Cairo.lineTo coord 1
		Cairo.stroke

		-- label grid lines
		spaceExtent <- Cairo.textExtents [' ']
		fe <- Cairo.fontExtents
		labels <- traverse sequence [(d, sizedText (dtbLabel dtbX d)) | d <- ds]

		let maxWidth = maximum $ map (Cairo.textExtentsXadvance . stExtents . snd) labels
		    xScale = 0.1/(Cairo.textExtentsXadvance spaceExtent + maxWidth)
		    yScale = 0.1/(Cairo.fontExtentsAscent fe + Cairo.fontExtentsDescent fe)
		    scale = min xScale yScale
		    halfScale = scale/2
		    y = 1.1 + (Cairo.fontExtentsAscent fe - Cairo.fontExtentsDescent fe) * halfScale

		scaleFontMatrix scale
		forM_ labels $ \(d, st) -> do
			let dx = Cairo.textExtentsXadvance (stExtents st) * halfScale
			Cairo.moveTo (xPos d - dx) y
			Cairo.showText (dtbLabel dtbX d)

	when hasState $ do
		-- draw grid lines
		Cairo.setDash [1/60, 3/60, 1/60, 1/60] (-1/120)
		flip M.traverseWithKey states $ \state ix -> when (isMajor state) $ do
			let y = 1 - ix / maxStateIx
			Cairo.moveTo 0 y
			Cairo.lineTo 1 y
		Cairo.stroke

		-- label grid lines
		valignText 1.05 0.2 . M.fromList $
			[ (yPos (State state), ("", "", state))
			| state <- M.keys states
			, isMajor state
			]

	when hasDiffTime $ do
		let ds = dtbGridLines dtbY
		-- draw grid lines
		Cairo.setDash [1/60, 1/60] (1/120)
		for_ ds $ \d -> do
			let coord = yPos (Time d)
			Cairo.moveTo 0 coord
			Cairo.lineTo 1 coord
		Cairo.stroke

		-- label grid lines
		valignText (-0.25) 0.2 . M.fromList $
			[(yPos (Time d), (dtbLabel dtbY d, "", "")) | d <- ds]

	when hasLogScale $ do
		-- draw grid lines
		Cairo.setDash [1/60, 5/60] (-5/120)
		for_ (lsbGridLines lsb) $ \n -> do
			let coord = yPos (LogScale n)
			Cairo.moveTo 0 coord
			Cairo.lineTo 1 coord
		Cairo.stroke

		-- label grid lines
		valignText (-0.4) 0.1 . M.fromList $
			[ (yPos (LogScale gridLine), logScaleShow gridLine)
			| gridLine <- lsbGridLines lsb
			]

	for_ (zip [0..] ptss) $ \(i, (lbl, pts)) -> do
		let dAngle = min (pi / numModules) (pi / 12)
		    angle = 2 * pi * i / numModules
		    RGB r g b = toSRGB (cieLAB d65 50 (100 * cos angle) (100 * sin angle))
		Cairo.setSourceRGB r g b
		for_ pts $ \pt -> do
			let xCoord = xPos (x pt)
			    yCoord = yPos (y pt)
			Cairo.moveTo xCoord yCoord
			Cairo.arc xCoord yCoord 0.03 (angle-dAngle) (angle+dAngle)
		Cairo.fill

valignText :: (Cairo.CairoString string, IsString string, Monoid string) =>
	Double -> Double -> Map Double (string, string, string) -> Render ()
valignText x w entries_ = do
	entries <- traverse valignInfo $
		if M.null entries_
		then M.singleton 0 ("", "", "")
		else entries_
	fe <- Cairo.fontExtents

	let maxPreWidth = maximum (vaiPrefixWidth <$> entries)
	    maxSufWidth = maximum (vaiSuffixWidth <$> entries)
	    maxWidth = maxPreWidth + maxSufWidth
	    xScale = w / maxWidth
	    yScale = chooseYScale xScale fe entries
	    scale = min xScale yScale
	    sepX = x + (w + (maxPreWidth - maxSufWidth)*scale)/2
	    dy = (Cairo.fontExtentsAscent fe - Cairo.fontExtentsDescent fe)*scale/2

	scaleFontMatrix scale

	void . flip M.traverseWithKey entries $ \y_ vai -> do
		let y = y_ + dy
		    sepWidth = scale * (Cairo.textExtentsXadvance . stExtents . vaiSeparator) vai
		Cairo.moveTo (sepX - vaiPrefixWidth vai*scale) y
		Cairo.showText (stString (vaiPrefix vai) <> stString (vaiSeparator vai))
		Cairo.moveTo (sepX + sepWidth/2) y
		Cairo.showText (stString (vaiSuffix vai))

scaleFontMatrix :: Double -> Render ()
scaleFontMatrix s = Cairo.setFontMatrix . CM.scalarMultiply s =<< Cairo.getFontMatrix

data SizedText string = SizedText
	{ stString :: string
	, stExtents :: Cairo.TextExtents
	}

sizedText :: Cairo.CairoString string => string -> Render (SizedText string)
sizedText s = SizedText s <$> Cairo.textExtents s

data VAlignInfo string = VAlignInfo
	{ vaiPrefix, vaiSeparator, vaiSuffix :: SizedText string
	, vaiPrefixWidth :: Double
	, vaiSuffixWidth :: Double
	}

valignInfo :: Cairo.CairoString string => (string, string, string) -> Render (VAlignInfo string)
valignInfo (pre_, sep_, suf_) = do
	pre <- sizedText pre_
	sep <- sizedText sep_
	suf <- sizedText suf_
	pure VAlignInfo
		{ vaiPrefix = pre
		, vaiSeparator = sep
		, vaiSuffix = suf
		, vaiPrefixWidth = dx pre + dx sep / 2
		, vaiSuffixWidth = dx suf + dx sep / 2
		}
	where dx = Cairo.textExtentsXadvance . stExtents

chooseYScale :: Double -> Cairo.FontExtents -> Map Double (VAlignInfo string) -> Double
chooseYScale def fe vai
	| M.size vai < 2 = def
	| otherwise = minDY / reqDY
	where
	minDY = minimum . ap (zipWith subtract) tail . M.keys $ vai
	reqDY = Cairo.fontExtentsHeight fe

-- TODO: the smallest (and largest?) state should always be major

data DiffTimeBounds = DiffTimeBounds
	{ dtbZero :: DiffTimeSpec
	, dtbGridLine :: DiffTimeSpec
	, dtbOne :: DiffTimeSpec
	, dtbLabel :: DiffTimeSpec -> String
	}

chooseDiffTimeBounds :: [DiffTimeSpec] -> DiffTimeBounds
chooseDiffTimeBounds [] = DiffTimeBounds
	{ dtbZero = 0
	, dtbGridLine = 1
	, dtbOne = 10
	, dtbLabel = \t -> show (floor t) ++ "s"
	}
chooseDiffTimeBounds ts = DiffTimeBounds
	{ dtbZero = smallAndRound
	, dtbGridLine = grid
	, dtbOne = bigAndRound
	, dtbLabel = label
	} where
	small = minimum ts
	big = maximum ts
	width = big - small
	gridTarget = width / 10
	(grid, label) = case M.lookupGE gridTarget roundIntervals of
		Just v -> v
		Nothing -> ( gridDays (60*60*24) (gridTarget / (60*60*24))
		           , \t -> let (pre, sep, suf) = logScaleShow . fromInteger . floor $ t/(60*60*24)
		                   in concat [pre, sep, suf, "d"]
		           )
	gridDays cur tgt
		| tgt > 5 = gridDays (cur*10) (tgt/10)
		| tgt > 2 = cur*5
		| tgt > 1 = cur*2
		| otherwise = cur
	smallAndRound = fromInteger (floor (small / grid)) * grid
	bigAndRound = fromInteger (ceiling (big / grid)) * grid

roundIntervals :: Map DiffTimeSpec (DiffTimeSpec -> String)
roundIntervals = M.fromList $ zip
	[      1,       2,       5,       10,       30 -- seconds
	,   60*1,    60*2,    60*5,    60*10,    60*30 -- minutes
	,60*60*1, 60*60*2, 60*60*6, 60*60*12           -- hours
	]
	([]
	++ replicate 5 (fmt "%m:%02S")
	++ replicate 5 (\t -> fmt (if t < 60*60 then "%-mm" else "%-hh%02Mm") t)
	++ replicate 4 (\t -> fmt (if t < 60*60*24 then "%-hh" else "%-dd%02Hh") t)
	)
	where
	fmt format = formatTime defaultTimeLocale format . ndt

	ndt :: DiffTimeSpec -> NominalDiffTime
	ndt = realToFrac

dtbGridLines :: DiffTimeBounds -> [DiffTimeSpec]
dtbGridLines dtb = [dtbZero dtb, dtbZero dtb + dtbGridLine dtb .. dtbOne dtb]

dtbPos :: DiffTimeBounds -> DiffTimeSpec -> Double
dtbPos dtb d = if o == z then 0.5 else realToFrac ((d - z) / (o - z)) where
	o = dtbOne dtb
	z = dtbZero dtb

data LogScaleBounds = LogScaleBounds
	{ lsbZero :: Double
	, lsbScale :: Double
	, lsbOffset :: Double
	, lsbGridLines :: [Double]
	} deriving (Eq, Ord, Read, Show)

chooseLogScaleBounds :: [Double] -> LogScaleBounds
chooseLogScaleBounds ns
	| null ns = LogScaleBounds
		{ lsbZero = 1
		, lsbScale = recip (log 10)
		, lsbOffset = 0
		, lsbGridLines = [1, 2, 5, 10]
		}
	| all (< 0) ns = let lsb = chooseLogScaleBounds (negate <$> ns) in LogScaleBounds
		{ lsbZero = -lsbZero lsb
		, lsbScale = -lsbScale lsb
		, lsbGridLines = negate <$> lsbGridLines lsb
		, lsbOffset = 1
		}
	| any (<= 0) ns = chooseLogScaleBounds (filter (>0) ns)
	| otherwise = LogScaleBounds
		{ lsbZero = zero
		, lsbScale = recip (log (one/zero))
		, lsbOffset = 0
		, lsbGridLines = gridLines
		} where
		small = minimum ns
		big = maximum ns

		oom x = floor (logBase 10 x)
		logMags@[smallLogMag, bigLogMag] = map oom [small, big]
		mags@[smallMag, bigMag] = map (10^^) logMags
		[smallRatio, bigRatio] = zipWith (/) [small, big] mags

		smallRoundRatio = fromMaybe 1 (S.lookupLE smallRatio roundDigits)
		bigRoundRatio = fromMaybe 10 (S.lookupGE bigRatio roundDigits)

		zero = smallRoundRatio * smallMag
		one = bigRoundRatio * bigMag

		allGrids = S.toAscList . S.fromList $
			[ gridLine
			-- bigLogMag is sometimes too small by 1 due to rounding, so
			-- include 10 as a choice of coefficient
			| ratio <- 10 : S.toList roundDigits
			, mag <- [smallLogMag..bigLogMag]
			, let gridLine = ratio*10^^mag
			, zero <= gridLine && gridLine <= one
			]
		gridLines = chooseAFewEvenly allGrids

chooseAFewEvenly :: [a] -> [a]
chooseAFewEvenly as
	| lines <= 6 = as
	| otherwise = head $ []
		++ [ every (gaps `quot` n) as
		   | n <- [5, 4, 3]
		   , gaps `rem` n == 0
		   ]
		++ [ indices [ (numerator * (lines-1) + 2) `quot` 4
		             | numerator <- [0 .. 4]
		             ]
		   ]
	where
	lines = length as
	gaps = lines-1

	every n [] = []
	every n (x:xs) = x : every n (drop (n-1) xs)

	indices = go 0 as where
		go _ [] _ = []
		go _ _ [] = []
		go i (x:xs) is@(i':is') = case compare i i' of
			LT -> go (i+1) xs is
			EQ -> x : go (i+1) xs is'
			-- should be impossible, but...
			GT -> if i' < 0 then go i (x:xs) is' else go 0 as is

roundDigits :: Set Double
roundDigits = S.fromList [1, 2, 5]

logScalePos :: LogScaleBounds -> Double -> Double
logScalePos lsb n = lsbOffset lsb + log (n/lsbZero lsb) * lsbScale lsb

logScaleShow :: Double -> (String, String, String)
logScaleShow = id
	. snd
	. head
	. sortOn (\(penalty, (pre, mid, suf)) -> (penalty + length (pre++mid++suf), mid))
	-- penalty: non-scientific notation is a bit nicer when it's possible, even
	-- to the point that we might want to choose it even when it's longer
	. (\v -> [(0, showDigitsFlat v), (1, showDigitsScientific v)])
	. dropZeros
	. sigFigs

sigFigs :: Double -> (Int, Int)
sigFigs 0 = (0, 0)
sigFigs n = (round (n/scale), magnitude) where
	-- 0.01 is a fudge factor for slightly-too-low log calculations
	magnitude = floor (logBase 10 (abs n) + 0.01) - 2
	scale = 10^^magnitude

dropZeros :: (Int, Int) -> (Int, Int)
dropZeros v@(mantissa, exponent) = case mantissa `quotRem` 10 of
	(0, 0) -> v
	(mantissa', 0) -> dropZeros (mantissa', exponent+1)
	_ -> v

showDigitsFlat :: (Int, Int) -> (String, String, String)
showDigitsFlat (mantissa, exponent)
	| exponent >= 0 = (show mantissa ++ replicate exponent '0', "", "")
	| otherwise = go mantissa exponent ""
	where
	go 0 e suf = (['-' | mantissa < 0] ++ "0", ".", replicate (-e) '0' ++ suf) -- fast path
	go m 0 suf = (show m, ".", suf)
	go m e suf = let (q, r) = m `quotRem` 10 in go q (e+1) (intToDigit (abs r) : suf)

showDigitsScientific :: (Int, Int) -> (String, String, String)
showDigitsScientific (mantissa, exponent) = go mantissa exponent "" where
	go m e suf
		| abs m < 10 = (show m ++ ['.' | not (null suf)] ++ suf, "e", show e)
		| otherwise = let (q, r) = m `quotRem` 10 in go q (e+1) (intToDigit (abs r) : suf)

-- The panedPosition offered by the library seems to be implemented in some
-- roundabout way that doesn't let notifyProperty work right.
panedPosition' :: PanedClass self => Attr self Int
panedPosition' = newAttrFromIntProperty "position"

updateTimerLabel :: Context -> IO ()
updateTimerLabel ctx = postGUIAsync $ do
	status <- readMVar (ctxUITimerLabelStatus ctx)
	case status of
		RunningSince from -> go runningAttrs . (`diffTimeSpec` from) =<< getCurrentTimeSpec
		Stopped dur -> go (if dur == 0 then commonAttrs else stoppedAttrs) dur
	where
	commonAttrs = tail [undefined
		, AttrWeight { paWeight = WeightHeavy, paStart = setLater, paEnd = setLater }
		, AttrSize { paSize = 32, paStart = setLater, paEnd = setLater }
		, AttrFamily { paFamily = "FreeMono,Latin Modern Mono,Noto Sans Mono CJK JP,DejaVu Sans Mono", paStart = setLater, paEnd = setLater }
		]
	runningAttrs = commonAttrs ++ tail [undefined
		, AttrForeground { paColor = Color 0x1000 0x8000 0x0000, paStart = setLater, paEnd = setLater }
		]
	stoppedAttrs = commonAttrs ++ tail [undefined
		, AttrForeground { paColor = Color 0xb000 0x2800 0x3800, paStart = setLater, paEnd = setLater }
		]
	setLater = error "updateTimerLabel attempted to use an attribute without setting its range first. This is a bug."

	go attrs dur = do
		mag <- readIORef (ctxTimeMagnitude ctx)
		let (mag', text) = showDiffTimeSpec dur mag
		labelSetText (ctxUITimerLabel ctx) text
		labelSetAttributes (ctxUITimerLabel ctx) [attr { paStart = 0, paEnd = length text } | attr <- attrs]
		when (mag /= mag') $ do
			writeIORef (ctxTimeMagnitude ctx) mag'
			redrawAllIntervals ctx

updateGraph :: Context -> IO ()
updateGraph = postGUIAsync . widgetQueueDraw . ctxUIGraph

-- 2^63 microseconds is about 300k years, so even though this technically can't
-- represent all the possible TimeSpec differences, it's gonna be good enough
-- for this application.

-- TODO: Read and Show instances that put the decimal point in the right place
newtype DiffTimeSpec = DiffTimeSpec { getMicros :: Int64 }
	deriving (Bounded, Eq, Ord, DB.FromField, DB.ToField, NFData)

-- | @fromInteger 1@ is one second
instance Num DiffTimeSpec where
	DiffTimeSpec a + DiffTimeSpec b = DiffTimeSpec (a+b)
	DiffTimeSpec a - DiffTimeSpec b = DiffTimeSpec (a-b)
	DiffTimeSpec a * DiffTimeSpec b = DiffTimeSpec (fromInteger ((toInteger a * toInteger b + 500000) `div` 1000000))
	negate (DiffTimeSpec a) = DiffTimeSpec (negate a)
	abs (DiffTimeSpec a) = DiffTimeSpec (abs a)
	signum (DiffTimeSpec a) = DiffTimeSpec (1000000 * signum a)
	fromInteger n = DiffTimeSpec (fromInteger (1000000 * n))

instance Real DiffTimeSpec where
	toRational (DiffTimeSpec a) = toRational a / 1000000

instance Fractional DiffTimeSpec where
	DiffTimeSpec a / DiffTimeSpec b = DiffTimeSpec (fromInteger ((a'*1000000 + b' `div` 2) `div` b')) where
		a' = toInteger a
		b' = toInteger b
	recip (DiffTimeSpec a) = DiffTimeSpec (fromInteger ((1000000000000 + a' `div` 2) `div` a')) where
		a' = toInteger a
	fromRational n = DiffTimeSpec (round (n*1000000))

instance RealFrac DiffTimeSpec where
	properFraction (DiffTimeSpec a) = (fromIntegral q, DiffTimeSpec r) where
		(q, r) = a `quotRem` 1000000

instance Enum DiffTimeSpec where
	succ (DiffTimeSpec a) = if a > maxBound - 1000000
		then error $ "DiffTimeSpec.succ: bad argument " ++ show a
		else DiffTimeSpec (a+1000000)
	pred (DiffTimeSpec a) = if a < minBound + 1000000
		then error $ "DiffTimeSpec.pred: bad argument " ++ show a
		else DiffTimeSpec (a-1000000)

	toEnum = coerce (toEnum :: Int -> Int64)
	fromEnum = coerce (int64 fromEnum)

	enumFrom (DiffTimeSpec a) = coerce [a, a+1000000 ..]
	enumFromThen = coerce (int64 enumFromThen)
	enumFromTo (DiffTimeSpec a) (DiffTimeSpec b) = coerce [a, a+1000000 .. b]
	enumFromThenTo = coerce (int64 enumFromThenTo)

instance Show DiffTimeSpec where
	show (DiffTimeSpec a) = concat . tail $ [undefined
		, if a < 0 then "-" else ""
		, show (abs q)
		, "."
		, pad (show (abs r))
		] where
		(q, r) = a `quotRem` 1000000
		pad s = replicate (6-length s) '0' ++ s

instance Read DiffTimeSpec where
	readPrec = readNumber $ \case
		Number n -> maybe (fail "not a number") convertDTS $ numberToFixed 6 n
		_ -> fail "not a number"
		where
		convertDTS (secs, micros) = pure . DiffTimeSpec . fromInteger $
			secs*1000000 + signum secs*micros

{-# INLINE int64 #-}
int64 :: (Int64 -> a) -> Int64 -> a
int64 f = f

diffTimeSpec :: TimeSpec -> TimeSpec -> DiffTimeSpec
diffTimeSpec (TimeSpec s ns) (TimeSpec s' ns') = DiffTimeSpec $
	1000000*(s-s') + (ns-ns'+500)`div`1000

addTimeSpec :: TimeSpec -> DiffTimeSpec -> TimeSpec
addTimeSpec ts (DiffTimeSpec micros) = ts + fromIntegral (1000*micros)

showDiffTimeSpec :: DiffTimeSpec -> TimeMagnitude -> (TimeMagnitude, String)
showDiffTimeSpec = showNominalDiffTime . realToFrac

showNominalDiffTime :: NominalDiffTime -> TimeMagnitude -> (TimeMagnitude, String)
showNominalDiffTime dur mag_ = (mag, formatTime defaultTimeLocale (formatString mag) dur) where
	mag = max mag_ (timeMagnitude dur)

timeMagnitude :: NominalDiffTime -> TimeMagnitude
timeMagnitude t
	| t < 60 = Seconds
	| t < 60*60 = Minutes
	| t < 24*60*60 = Hours
	| otherwise = Days (go (10*24*60*60) 1)
	where
	go t' d
		| t < t' = d
		-- technically could overflow, but like... 10^255 days is a really lot of days
		| otherwise = go (10*t') (1+d)

formatString :: TimeMagnitude -> String
formatString Seconds = "%03ES"
formatString Minutes = "%02M:" ++ formatString Seconds
formatString Hours = "%02H:" ++ formatString Minutes
formatString (Days n) = "%0" ++ show n ++ "dd " ++ formatString Hours

-- TODO
redrawAllIntervals :: Context -> IO ()
redrawAllIntervals _ = pure ()

getCurrentTimeSpec :: IO TimeSpec
getCurrentTimeSpec = getTime MonotonicRaw

stdinThread :: Context -> IO loop
stdinThread (ctxInput -> chan) = forever $ do
	line <- T.getLine
	now <- getCurrentTimeSpec
	writeChan chan (now, line)

data Event = Event
	{ eFrame :: Maybe Integer
	, eMicrosecond :: DiffTimeSpec
	, eState :: Text
	} deriving (Eq, Ord, Read, Show)

data Command = StateChange Event | Stop DiffTimeSpec deriving (Eq, Ord, Read, Show)

stop :: Text
stop = "STOP"

-- TODO: catch EOF and do something sensible (probably just quit?)
-- TODO: catch DB connection errors
parserThread :: Context -> IO loop
parserThread ctx = DB.connectPostgreSQL (db (ctxConfig ctx)) >>= go where
	go conn = case pFPS (ctxProfile ctx) of
		Just{} -> idling
		Nothing -> noFramesRunning Nothing Nothing
		where
		getEvent = do
			(now, t) <- readChan (ctxInput ctx)
			let (frame_, preState) = T.break (' '==) t
			    state = T.drop 1 preState
			case (T.signed T.decimal frame_, T.null preState) of
				(Right (frame, T.null -> True), False)
					| state == stop -> pure $ Left (Just now)
					| ctxValidState ctx state -> pure $ Right (frame, now, state)
					| otherwise -> Left Nothing <$ writeChan (ctxErrors ctx)
						(printf "Ignoring unknown state %s at frame %d, time %s" (show state) frame (show now))
				_ -> Left Nothing <$ writeChan (ctxErrors ctx) ("Ignoring malformed state change " ++ show t)

		idling = getEvent >>= \case
			Right (frame, now, state) -> do
				broadcast now Event { eFrame = Just 0, eMicrosecond = 0, eState = state }
					Nothing
					(running (negate frame) frame)
					stopping
			Left (Just now) -> warnIgnoreStop >> idling
			Left Nothing -> idling

		running dFrame frame runStart splitKey = getEvent >>= \case
			Right (frame', now, state)
				| frame' >= frame -> do
					broadcast runStart
						Event
							{ eFrame = Just (dFrame+frame')
							, eMicrosecond = diffTimeSpec now runStart
							, eState = state
							}
						(Just splitKey)
						(running dFrame frame')
						stopping
				| otherwise -> do
					let ddFrame = roundUpToPowerOf2 (max frame (if signed then -1-frame' else frame))
					    signed = frame' < 0
					    dFrame' = dFrame + if signed then shiftL ddFrame 1 else ddFrame
					writeChan (ctxErrors ctx) $ printf
						"New frame %d is smaller than old frame %d. Assuming this is due to %s overflow at %d."
						frame'
						frame
						(if signed then string "signed" else string "unsigned")
						ddFrame
					broadcast runStart
						Event
							{ eFrame = Just (dFrame'+frame')
							, eMicrosecond = diffTimeSpec now runStart
							, eState = state
							}
						(Just splitKey)
						(running dFrame' frame')
						stopping
			Left (Just now) -> stopTimer (diffTimeSpec now runStart) >> idling
			Left Nothing -> running dFrame frame runStart splitKey

		stopping = getEvent >>= \case
			Left (Just now) -> idling
			_ -> stopping

		-- invariant: arguments are both Just or both Nothing
		noFramesRunning runStartM splitKey = do
			(now, state) <- readChan (ctxInput ctx)
			let runStart = fromMaybe now runStartM
			case (state == stop, ctxValidState ctx state) of
				(True, _) -> do
					maybe warnIgnoreStop (stopTimer . diffTimeSpec now) runStartM
					noFramesRunning Nothing Nothing
				(_, True) -> broadcast runStart
					Event
						{ eFrame = Nothing
						, eMicrosecond = diffTimeSpec now runStart
						, eState = state
						}
					splitKey
					(\rs sk -> noFramesRunning (Just rs) (Just sk))
					noFramesStopping
				_ -> do
					writeChan (ctxErrors ctx)
						(printf "Ignoring unknown state %s at time %s" (show state) (show now))
					noFramesRunning runStartM splitKey

		noFramesStopping = do
			(now, state) <- readChan (ctxInput ctx)
			if state == stop
			then noFramesRunning Nothing Nothing
			else noFramesStopping

		broadcast runStart event splitKey continue won = do
			-- tell the database
			(runID, seqNo) <- case splitKey of
				Just (runID, seqNo) -> pure (runID, seqNo+1)
				Nothing -> do
					[DB.Only runID] <- DB.query conn
						"insert into run (game, fps, started) values (?, ?, now()) returning id"
						(pGame (ctxProfile ctx), pFPS (ctxProfile ctx))
					pure (runID :: ID, 0 :: ID)
			DB.execute conn
				"insert into split (run, seq_no, state, microsecond, frame) values (?, ?, ?, ?, ?)"
				(runID, seqNo, eState event, eMicrosecond event, eFrame event)

			-- tell the UI
			atomically $ do
				sr <- readTVar (ctxStateRange ctx)
				case srInsert (pSortOrder (ctxProfile ctx)) (eState event) sr of
					(Any True, sr') -> writeTVar (ctxStateRange ctx) sr'
					_ -> pure ()
			status <- takeMVar (ctxUITimerLabelStatus ctx)
			putMVar (ctxUITimerLabelStatus ctx) $ case status of
				RunningSince{} -> status
				_ -> RunningSince runStart
			postGUIAsync $ do
				case status of
					RunningSince{} -> pure ()
					_ -> treeStoreClear (ctxEventStore ctx)
				logEvent ctx event

			-- tell the modules
			writeChan (ctxModuleInput ctx) (StateChange event)

			-- tell the parser thread
			if eState event == pTarget (ctxProfile ctx)
			then stopTimer (eMicrosecond event) >> won
			else continue runStart (runID, seqNo)

		stopTimer dur = do
			-- nothing interesting to tell the database
			-- tell the UI
			status <- takeMVar (ctxUITimerLabelStatus ctx)
			putMVar (ctxUITimerLabelStatus ctx) $ case status of
				RunningSince{} -> Stopped dur
				_ -> status
			-- tell the modules
			writeChan (ctxModuleInput ctx) (Stop dur)
			-- the parser thread already knows

		warnIgnoreStop = writeChan (ctxErrors ctx) "Ignoring redundant STOP"

ctxValidState :: Context -> Text -> Bool
ctxValidState ctx s = maybe True (O.member s) (ctxAllStates ctx)

roundUpToPowerOf2 :: Integer -> Integer
roundUpToPowerOf2 n
	| n <= 1 = 1
	| otherwise = go 1 2
	where
	go exp pow
		| n <= pow = head [bit (e+1) | e <- [exp, exp-1 .. exp `rem` 2], testBit n e]
		| otherwise = go (2*exp) (pow*pow)

logEvent :: Context -> Event -> IO ()
logEvent ctx e = do
	n <- treeModelIterNChildren store Nothing
	case pMajorStates (ctxProfile ctx) of
		Just major | eState e `S.member` major -> do
			es <- removeMinorEvents major [] (n-1)
			n' <- treeModelIterNChildren store Nothing
			treeStoreInsertTree store [] n' $ Node e [Node e' [] | e' <- es]
		_ -> treeStoreInsert store [] n e
	where
	store = ctxEventStore ctx
	removeMinorEvents major minor n
		| n < 0 = pure minor
		| otherwise = do
			e' <- treeStoreGetValue store [n]
			if eState e' `S.member` major
			then pure minor
			else do
				treeStoreRemove store [n]
				removeMinorEvents major (e':minor) (n-1)

class Ord a => IsText a where
	fromText :: Text -> a
	toText :: a -> Text

instance IsText Text where
	fromText = id
	toText = id

instance IsText HumanOrder where
	fromText = humanOrder
	toText = unHumanOrder

instance IsText a => IsText (Down a) where
	fromText = Down . fromText
	toText (Down a) = toText a

data StateCache where
	Ordered :: OSet Text -> StateCache
	Sorted :: IsText a => Set a -> StateCache
	Uninitialized :: StateCache

stateCache :: Context -> [Text] -> StateCache
stateCache ctx ss = case pSortOrder (ctxProfile ctx) of
	Ascending -> Sorted . S.fromList . map humanOrder $ ss
	Descending -> Sorted . S.fromList . map (Down . humanOrder) $ ss
	AscendingOn o -> Ordered o

data StateRange = StateRange
	{ srStateCache :: StateCache
	, srMin :: HumanOrder
	, srMax :: HumanOrder
	}

srMinT, srMaxT :: StateRange -> Text
srMinT = toText . srMin
srMaxT = toText . srMax

srInsert :: SortOrder -> Text -> StateRange -> (Any, StateRange)
srInsert so t_ sr = pure StateRange
	<*> case srStateCache sr of
		Sorted ts -> (Any (S.notMember t' ts), Sorted (S.insert t' ts)) where t' = fromText t_
		_ -> (Any False, srStateCache sr) -- slightly optimistic... if it's Ordered, assumes t appears in the OSet
	<*> (Any (soLt so t (srMin sr)), soMin so t (srMin sr))
	<*> (Any (soLt so (srMax sr) t), soMax so t (srMax sr))
	where t = humanOrder t_

ctxStates :: Context -> IO [Text]
ctxStates ctx = do
	sr <- readTVarIO (ctxStateRange ctx)
	case srStateCache sr of
		Ordered ss -> case (O.findIndex (srMinT sr) ss, O.findIndex (srMaxT sr) ss) of
			(Just iMin, Just iMax) -> pure [fromJust (ss `O.elemAt` i) | i <- [iMin .. iMax]]
			_ -> failedRange sr $ printf "StateRange bounds %s-%s not found in StateCache %s"
				(show (srMinT sr))
				(show (srMaxT sr))
				(show ss)
		Sorted ss -> case (S.lookupIndex eMin ss, S.lookupIndex eMax ss) of
			(Just iMin, Just iMax) -> pure
				. map toText
				. S.toAscList
				. S.take (iMax-iMin+1)
				. S.drop iMin
				$ ss
			_ -> failedRange sr $ printf "StateRange bounds %s-%s not found in StateCache %s"
				(show (srMinT sr))
				(show (srMaxT sr))
				(show (S.map toText ss))
			where
			eMin = fromText (srMinT sr)
			eMax = fromText (srMaxT sr)
		Uninitialized -> do
			mCache <- initializeStateCache ctx
			cache <- case mCache of
				Nothing -> stateCache ctx <$>
					failedRange sr "Could not contact database to initialize state cache"
				Just cache -> pure cache
			atomically $ do
				sr' <- readTVar (ctxStateRange ctx)
				writeTVar (ctxStateRange ctx) (sr' { srStateCache = cache })
			ctxStates ctx
	where
	failedRange sr err = do
		writeChan (ctxErrors ctx) err
		pure $ [srMinT sr | srMin sr /= srMax sr] ++ [srMaxT sr]

initializeStateCache :: Context -> IO (Maybe StateCache)
initializeStateCache ctx = case pSortOrder (ctxProfile ctx) of
	AscendingOn ss -> pure (Just (Ordered ss))
	Ascending -> fetchFromDB
	Descending -> fetchFromDB
	where
	fetchFromDB = do
		mConn <- try (DB.connectPostgreSQL (db (ctxConfig ctx)))
		case mConn of
			Left err -> pure (const Nothing (DB.sqlState err {- type-checking hint -}))
			Right conn -> do
				ss <- DB.query conn
					"select distinct state from run join split on run.id = split.run where run.game = ?"
					(DB.Only (pGame (ctxProfile ctx)))
				DB.close conn
				pure . Just . stateCache ctx . fromOnlys $ ss

fromOnlys :: [DB.Only a] -> [a]
fromOnlys = coerce

heartbeatThread :: Context -> IO loop
heartbeatThread ctx = forever $ do
	updateTimerLabel ctx
	threadDelay (1000000`div`30) -- update at 30fps...ish

data Dependent
	= State Text
	| Time DiffTimeSpec
	| LogScale Double
	deriving (Eq, Ord, Read, Show)

instance NFData TimeSpec where
	rnf (TimeSpec secs nsecs) = rnf secs `seq` rnf nsecs

instance NFData Dependent where
	rnf (State s) = rnf s
	rnf (Time t) = rnf t
	rnf (LogScale n) = rnf n

data ModulePoint = ModulePoint
	{ x :: DiffTimeSpec
	, y :: Dependent
	} deriving (Eq, Ord, Read, Show)

instance NFData ModulePoint where
	rnf (ModulePoint xVal yVal) = rnf xVal `seq` rnf yVal

-- | Common case: we've just computed an offset to add to an event's time.
dtDependent :: Event -> DiffTimeSpec -> Dependent
dtDependent e dt = Time (eMicrosecond e + dt)

-- | Common case: we've just computed an offset to add to an event's time.
dtPoint :: Event -> DiffTimeSpec -> ModulePoint
dtPoint Event { eMicrosecond = micro } dt = ModulePoint micro (Time (micro + dt))

data Module = Module
	{ mLabel :: String
	, mLaunch :: Context -> Chan Command -> TVar [ModulePoint] -> IO Void
	}

-- | If your module only ever adds points, and only needs to see the latest
-- event to know which points it wants to add, you can use this wrapper to do
-- some of the bookkeeping for you.
appendOnlyModule :: String -> (Context -> IO a) -> (a -> Command -> IO [ModulePoint]) -> Module
appendOnlyModule label initialize newPoints = Module label $ \ctx i o -> do
	a <- initialize ctx
	let go pts = do
	    	me <- readChan i
	    	dpts <- newPoints a me
	    	let pts' = dpts ++ pts
	    	when (not (null dpts)) $ do
	    		atomically (writeTVar o pts')
	    		updateGraph ctx
	    	go $ case me of
	    		Stop{} -> []
	    		StateChange{} -> pts'
	go []

-- | If you meet the requirements of 'appendOnlyModule', and additionally your
-- new points always occur at the same time as the most recent event, you can
-- use this to add the x-coordinate for you.
synchronousOnlyModule :: String -> (Context -> IO a) -> (a -> Command -> IO [Dependent]) -> Module
synchronousOnlyModule label initialize newDependents
	= appendOnlyModule label initialize $ \a me ->
		map (ModulePoint (meMicrosecond me)) <$> newDependents a me
	where
	meMicrosecond (StateChange e) = eMicrosecond e
	meMicrosecond (Stop t) = t

stateProgressModule :: Module
stateProgressModule = synchronousOnlyModule "state progress" def
	$ \ ~() me -> pure [State (eState e) | StateChange e <- [me]]

pbModule :: Module
pbModule = synchronousOnlyModule "PB" initialize addPBRun
	where
	initialize ctx = do
		pbRef <- newIORef Nothing
		eConn <- try (DB.connectPostgreSQL (db (ctxConfig ctx)))
		case eConn of
			Left err -> writeChan (ctxErrors ctx) $ printf
				"Error connecting to database in pbModule: %s"
				(displayException @DB.SqlError err)
			Right conn -> pure ()
		pure (pbRef, eConn, pGame (ctxProfile ctx), pTarget (ctxProfile ctx))

	addPBRun (_, Left err, _, _) _ = pure []
	addPBRun (pbRef, _, _, _) Stop{} = [] <$ writeIORef pbRef Nothing
	addPBRun a@(pbRef, Right conn, game, tgt) me@(StateChange e) = do
		mPB <- readIORef pbRef
		case mPB of
			Just pb -> pure $ case M.lookup (eState e) pb of
				Nothing -> []
				Just t -> [dtDependent e ((pb M.! tgt) - t)]
			Nothing -> do
				runs <- DB.query @_ @(ID, DiffTimeSpec) conn findPBRunQuery (eState e, tgt, game)
				pb <- case runs of
					[] -> pure M.empty
					(runID, _):_ -> M.fromList <$> DB.query conn
						"select state, min(microsecond) from split where run = '?' group by state"
						(DB.Only runID)
				writeIORef pbRef (Just pb)
				addPBRun a me

	findPBRunQuery = "\
		\select b_run, e_microsecond - b_microsecond as duration \
		\from \
			\(\
				\(\
					\select run as b_run, microsecond as b_microsecond, seq_no as b_seq_no \
					\from split \
					\where state = ? \
				\) as beginning \
				\join \
				\( \
					\select run as e_run, microsecond as e_microsecond, seq_no as e_seq_no \
					\from split \
					\where state = ? \
				\) as ending \
				\on b_run = e_run and b_seq_no < e_seq_no \
			\) as interval \
			\join \
			\run \
			\on interval.b_run = run.id \
		\where run.game = ? \
		\order by duration, b_run \
		\limit 1"

data Split = Split
	{ run :: ID
	, seqNo :: Int32
	, state :: Text
	, microsecond :: DiffTimeSpec
	, frame :: Maybe Int32
	} deriving (Eq, Ord, Read, Show)
	-- BEWARE: some functions assume that the Ord instance compares on run
	-- first, then seqNo, then other stuff

instance DB.FromRow Split where
	fromRow = pure Split
		<*> DB.field
		<*> DB.field
		<*> DB.field
		<*> DB.field
		<*> DB.field

allSplits :: Context -> Connection -> IO [Split]
allSplits (ctxProfile -> p) conn = do
	ids <- DB.query conn "select id from run where game = ?" (DB.Only (pGame p))
	DB.query conn "select run, seq_no, state, microsecond, frame from split where run in ?"
		. DB.Only
		. DB.In
		. fromOnlys @ID
		$ ids

data Segment = Segment
	{ endState :: Text
	, duration :: DiffTimeSpec
	, frames :: Maybe Int32
	} deriving (Eq, Ord, Read, Show)

finishableSegments :: Context -> [Split] -> Map Text [Segment]
finishableSegments ctx splits = id
	. M.fromListWith (++)
	. concat
	. ap (zipWith mkSegment) tail
	$ sorted
	where
	sorted = sort splits
	finishable = unsafeFinishableStates ctx sorted
	mkSegment split split' =
		[(state split, [Segment
			{ endState = state split'
			, duration = microsecond split' - microsecond split
			, frames = liftA2 (-) (frame split') (frame split)
			}])
		| run split == run split'
		, seqNo split + 1 == seqNo split'
		, state split `S.member` finishable
		, state split' `S.member` finishable
		]

-- unsafe because it assumes the provided split list is already sorted
unsafeFinishableStates :: Context -> [Split] -> Set Text
unsafeFinishableStates (ctxProfile -> p) = id
	. restrict
	. dfs (pTarget p)
	. M.fromListWith (++)
	. concat
	. ap (zipWith predecessor) tail
	where
	-- paranoia: what if there are two profiles with the same game name? (but
	-- not quite paranoid enough to deal with two profiles with the same game
	-- name and shared states. at some point you finally have to give in to
	-- GIGO)
	restrict ss = case soAllStates (pSortOrder p) of
		Just ss' -> S.intersection (O.toSet ss') ss
		Nothing -> ss

	predecessor s s' =
		[(state s', [state s]) | run s == run s' && seqNo s + 1 == seqNo s']

dfs :: Ord a => a -> Map a [a] -> Set a
dfs root neighbors = go root (S.singleton root) where
	go s seen = case M.lookup s neighbors of
		Nothing -> seen
		Just (S.fromList -> ss) -> foldr go (ss `S.union` seen) (ss S.\\ seen)

-- This produces the expected (in the statistics sense of the term) run time
-- from all source states that can reach the target to the target. The model
-- used here is that if you are currently in state S0, then the next state you
-- arrive at S1 and the time it takes to get there t are sampled IID from some
-- (unknown) distribution.
--
-- Suppose for a moment that we *did* know that distribution for each state. Then
--
-- E[T → T] = 0
--
-- E⎡S  → T⎤ = ∑ p  * ⎛t   + E[S  → T]⎞  	whenever S  ≠ T
--  ⎣ i    ⎦   j  ij  ⎝ ij      j     ⎠  	          i
--
-- where T is the target state, S  ranges over states that can reach T (as does
--                               i
--
-- S ⎞, p   is the probability of transitioning to state S  given that you are
--  j⎠   ij                                               j
--
-- in state S  right now, and t   is the expected duration of a transition from
--           i                 ij
--
-- S  to S .
--  i     j
--
-- Now in the case where we do *not* know the distribution for each state, we
-- can approximate it by the uniform distribution over transition/duration
-- pairs from previous runs. This gives us estimates of p   and t   and a
--                                                       ij      ij
-- collection of linear equations with E⎡S  → T⎤ as the variables. We can
--                                      ⎣ i    ⎦
-- estimate these expectations by solving the system of linear equations.
expectedRunTimes :: Context -> Connection -> IO (Map Text DiffTimeSpec)
expectedRunTimes ctx conn = do
	splits <- allSplits ctx conn
	let target = pTarget (ctxProfile ctx)
	    stats = summarize <$> finishableSegments ctx splits
	    ss = S.insert target (M.keysSet stats)
	    nonTargets = S.delete target ss -- N.B. (S.delete k . S.insert k) is not necessarily id
	    -- this does all the addition in the world of exact integers, then at
	    -- the last second converts to Double; this avoids rounding issues to
	    -- the extent possible (at the unlikely cost of overflow issues)
	    summarize segs = (\(n, dur) -> (n, realToFrac dur)) <$> M.fromListWith pointwiseAddition
	    	[(endState seg, (1, duration seg)) | seg <- segs]
	    pointwiseAddition (n, dur) (n', dur') = (n+n', dur+dur')
	    outEdgeCounts = sum . fmap fst <$> stats
	    nonTargetsList = S.toList nonTargets
	    transitionMatrix = Matrix.matrix (S.size nonTargets)
	    	[ M.findWithDefault 0 to outEdges / outEdgeCount
	    	| from <- nonTargetsList
	    	, let outEdges = fst <$> M.findWithDefault M.empty from stats
	    	      outEdgeCount = max 1 (M.findWithDefault 1 from outEdgeCounts)
	    	, to <- nonTargetsList
	    	]
	    edgeTimes = Matrix.vector
	    	[ sum
	    		[ M.findWithDefault 0 to outEdges / outEdgeCount
	    		| to <- S.toList ss
	    		]
	    	| from <- nonTargetsList
	    	, let outEdges = snd <$> M.findWithDefault M.empty from stats
	    	      outEdgeCount = max 1 (M.findWithDefault 1 from outEdgeCounts)
	    	]
	when (M.keysSet stats /= nonTargets) . writeChan (ctxErrors ctx) $ printf
		"Some finishable states didn't have any edges to other finishable states or something? Weird. (nonTargets, stats) = %s"
		(show (nonTargets, stats))
	let expectedTimes = (Matrix.ident (S.size nonTargets) - transitionMatrix) Matrix.<\> edgeTimes
	pure
		. M.insert target 0
		. M.fromList
		. zipWith (\s t -> (s, realToFrac t)) nonTargetsList
		$ Matrix.toList expectedTimes

expectedModule :: Module
expectedModule = synchronousOnlyModule "expected" initialize handleEvent where
	initialize ctx = do
		conn <- DB.connectPostgreSQL (db (ctxConfig ctx))
		expectedsRef <- newIORef Nothing
		pure (ctx, conn, expectedsRef)

	handleEvent (ctx, conn, expectedsRef) e = case e of
		Stop{} -> [] <$ writeIORef expectedsRef Nothing
		StateChange e -> do
			expectedsMaybe <- readIORef expectedsRef
			expecteds <- case expectedsMaybe of
				Just expecteds -> pure expecteds
				Nothing -> do
					expecteds <- expectedRunTimes ctx conn
					writeIORef expectedsRef (Just expecteds)
					pure expecteds
			pure . map (dtDependent e) . toList
				$ M.lookup (eState e) expecteds

type Runs = Map Text (Multiset DiffTimeSpec)

runGenerationThread :: Context -> Chan Command -> IO loop
runGenerationThread ctx cmds = do
	conn <- DB.connectPostgreSQL (db (ctxConfig ctx))
	waiting conn
	where
	runTVar = ctxRandomRuns ctx
	target = pTarget (ctxProfile ctx)
	maxRunCount = 5000 -- TODO: configurable

	waiting conn = do
		cmd <- readChan cmds
		case cmd of
			Stop{} -> waiting conn
			StateChange e -> do
				atomically $ modifyTVar runTVar (\(epoch, _) -> (epoch+1, M.empty))
				segs <- finishableSegments ctx <$> allSplits ctx conn
				generating conn segs (eState e)

	generating conn segs state = do
		(_, runs) <- readTVarIO runTVar
		let remainingRunCount = maxRunCount - if state `M.member` segs
		    	then fromIntegral (MS.size (M.findWithDefault MS.empty state runs))
		    	else maxRunCount
		replicateM_ remainingRunCount $ do
			events <- randomRun segs target state
			case events of
				(s,_):_ | s == target -> atomically . modifyTVar' runTVar $
					\(epoch, runs) -> (,) (epoch+1) $! addRun events runs
				_ -> pure ()
		cmd <- readChan cmds
		case cmd of
			Stop{} -> waiting conn
			StateChange e -> generating conn segs (eState e)

randomRun :: Map Text [Segment] -> Text -> Text -> IO [(Text, DiffTimeSpec)]
randomRun segMap_ tgt = go 0 [] where
	segMap = V.fromList <$> segMap_
	go t prev s
		| s == tgt = pure prev'
		| otherwise = case M.lookup s segMap of
			Just segs -> do
				seg <- (segs V.!) <$> randomRIO (0, V.length segs - 1)
				go (t+duration seg) prev' (endState seg)
			Nothing -> pure prev'
		where prev' = (s,t):prev

addRun :: [(Text, DiffTimeSpec)] -> Runs -> Runs
addRun [] runs = runs
addRun es@((_, t):_) runs = M.S.unionWith MS.union newRuns runs where
	-- const keeps the last (earliest) one in the list if we hit a given state
	-- multiple times. If we were paranoid, we could use max instead. We don't
	-- want to keep them both because they are correlated -- and so not IID
	-- samples.
	newRuns = M.fromListWith const
		[(state, MS.singleton (t-t')) | (state, t') <- es]

-- | Simple modules that want access to random runs generated by the
-- 'runGenerationThread' can use this to do the bulk of the bookkeeping for
-- them. They will automatically get incremental updating if generation takes a
-- long time.
withRandomRuns :: String ->
	-- | Initialize some custom data on module launch
	(Context -> IO a) ->
	-- | Initialize some custom data on the first event from a run
	(Context -> a -> Event -> IO b) ->
	-- | Generate some points, given the most up-to-date random runs
	(Context -> a -> b -> Runs -> Event -> IO [ModulePoint]) ->
	Module
withRandomRuns s initializeLaunch initializeRun generate = Module s $ \ctx i_ o -> do
	(epoch, _) <- readTVarIO (ctxRandomRuns ctx)

	i <- newChan
	forkIO . forever $ readChan i_ >>= writeChan i . Right
	forkIO $ awaitEpochChanges (ctxRandomRuns ctx) i epoch

	a <- initializeLaunch ctx
	waiting ctx i o a
	where
	awaitEpochChanges i_ i epoch = do
		(epoch', runs) <- atomically $ do
			v@(epoch', runs) <- readTVar i_
			when (epoch == epoch') retry
			pure v
		writeChan i (Left runs)
		awaitEpochChanges i_ i epoch'

	waiting ctx i o a = do
		modification <- readChan i
		case modification of
			Left{} -> waiting ctx i o a
			Right Stop{} -> waiting ctx i o a
			Right (StateChange e) -> do
				b <- initializeRun ctx a e
				running ctx i o a b M.empty [e]

	running ctx i o a b runs events = do
		modification <- readChan i
		case modification of
			Left runs' -> redraw ctx i o a b runs' events
			Right (StateChange event) -> redraw ctx i o a b runs (event:events)
			Right Stop{} -> waiting ctx i o a

	redraw ctx i o a b runs events = do
		pts <- concat <$> traverse (generate ctx a b runs) events
		pts `deepseq` atomically (writeTVar o pts)
		updateGraph ctx
		running ctx i o a b runs events

percentileModule :: Double -> Module
percentileModule p_ = withRandomRuns (pString ++ " percentile")
	(\_ -> pure ())
	(\_ _ _ -> pure ())
	$ \_ _ _ runs event -> pure . toList $ do
		dts <- M.lookup (eState event) runs
		dt <- dts MS.!? round ((1-p) * fromIntegral (MS.size dts - 1))
		pure $ dtPoint event dt
	where
	p = min 1 (max 0 p_)

	pString = if fromIntegral pBigInt == pBig
		then show pBigInt ++ ordinalSuffix pBigInt
		else show pBig ++ "th"
		where
		pBig = 100*p
		pBigInt = round pBig
	ordinalSuffix 1 = "st"
	ordinalSuffix 2 = "nd"
	ordinalSuffix 3 = "rd"
	ordinalSuffix _ = "th"

pbChanceModule :: Module
pbChanceModule = withRandomRuns "attempts to PB"
	(DB.connectPostgreSQL . db . ctxConfig)
	getPBInfo
	(\_ _ -> computeAttempts)
	where
	getPBInfo (ctxProfile -> p) conn e = do
		pbLengths <- DB.query conn
			"select ending.microsecond - beginning.microsecond as duration \
			\from (select id from run where game = ?) as this_game \
			\join (select run, seq_no, microsecond from split where state = ?) as beginning \
			\on this_game.id = beginning.run \
			\join (select run, seq_no, microsecond from split where state = ?) as ending \
			\on this_game.id = ending.run and beginning.seq_no <= ending.seq_no \
			\order by duration \
			\limit 1"
			(pGame p, eState e, pTarget p)
		pure (fromOnlys pbLengths)

	computeAttempts pbLengths runs event = pure $ do
		pbLength <- pbLengths
		Just dts <- [M.lookup (eState event) runs]
		let nLT = MS.countLT (pbLength - eMicrosecond event) dts
		guard (eMicrosecond event < pbLength)
		guard (nLT > 0)
		pure ModulePoint
			{ x = eMicrosecond event
			, y = LogScale $ fromIntegral (MS.size dts) / fromIntegral nLT
			}

-- TODO: make this configurable
allModules :: [Module]
allModules = [stateProgressModule, percentileModule 0.1, percentileModule 0.9, expectedModule, pbModule, pbChanceModule]

moduleThread :: Context -> IO loop
moduleThread ctx = do
	(chans, outs) <- fmap unzip . for allModules $ \m -> do
		chan <- newChan
		tvar <- newTVarIO []
		forkIO (() <$ mLaunch m ctx chan tvar)
		pure (chan, (mLabel m, tvar))
	atomically $ writeTVar (ctxModuleOutputs ctx) outs
	-- TODO: when modules are configurable, only launch this thread when it's
	-- needed
	chan <- newChan
	forkIO (runGenerationThread ctx chan)
	forever $ do
		mEvent <- readChan (ctxModuleInput ctx)
		for_ (chan:chans) $ \chan -> writeChan chan mEvent

-- TODO: should we have an MVar or something so the main thread can wait till
-- this finishes?
preferencesThread :: Context -> IO loop
preferencesThread ctx = do
	initialPreferences <- loadPreferences (argProfile (ctxArguments ctx))
	go $ fromMaybe impossiblePrefs initialPreferences
	where
	go prefs = do
		prefs' <- debounceReadTChanIO (ctxPreferences ctx)
		when (prefs /= prefs') (savePreferences ctx prefs')
		go prefs'

	impossiblePrefs = Preferences
		{ windowRectangle = Rectangle (-1) (-1) (-1) (-1)
		, splitHeight = -1
		, graphWidth = -1
		}

debounceReadTChanIO :: TChan a -> IO a
debounceReadTChanIO chan = atomically (readTChan chan) >>= go where
	go a = do
		threadDelay 5000000 -- 5-10 seconds after the last mouse movement ought to be soon enough
		(continue, a') <- atomically (emptyTChan chan a)
		(if continue then go else pure) a'

emptyTChan :: TChan a -> a -> STM (Bool, a)
emptyTChan chan = go False where
	go haveRead res = tryReadTChan chan >>= \case
		Nothing -> pure (haveRead, res)
		Just res' -> go True res'

savePreferences :: Context -> Preferences -> IO ()
savePreferences ctx prefs = do
	dir <- preferencesDir
	let fp = dir </> argProfile (ctxArguments ctx)
	    write = do
	    	createDirectoryIfMissing True dir
	    	writeFile fp (preferencesToString prefs)

	    complain :: IOException -> IO ()
	    complain err = do
	    	writeChan (ctxErrors ctx) $ printf
	    		"Error while attempting to write preferences to %s (will retry at next change): %s"
	    		fp (displayException err)
	catch write complain

preferencesToString :: Preferences -> String
preferencesToString = flip shows "\n" . C.pretty . preferencesToValue where
	preferencesToValue (Preferences (Rectangle x y w h) sh gw) = C.Sections ()
		[ C.Section () "window" $ C.Sections ()
			[ intSection "x" x
			, intSection "y" y
			, intSection "width" w
			, intSection "height" h
			]
		, C.Section () "splits" $ C.Sections () [intSection "height" sh]
		, C.Section () "graph" $ C.Sections () [intSection "width" gw]
		]
	intSection nm = C.Section () nm . C.Number () . C.integerToNumber . toInteger

max_stars_repo_name: dmwit/vertumnion

max_issues_repo_name: dmwit/vertumnion

max_forks_repo_name: dmwit/vertumnion


text:
{-# LANGUAGE FlexibleContexts #-}
-- |
-- Module    : System.Random.MWC.CondensedTable
-- Copyright : (c) 2012 Aleksey Khudyakov
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Table-driven generation of random variates.  This approach can
-- generate random variates in /O(1)/ time for the supported
-- distributions, at a modest cost in initialization time.
module System.Random.MWC.CondensedTable (
    -- * Condensed tables
    CondensedTable
  , CondensedTableV
  , CondensedTableU
  , genFromTable
    -- * Constructors for tables
  , tableFromProbabilities
  , tableFromWeights
  , tableFromIntWeights
    -- ** Disrete distributions
  , tablePoisson
  , tableBinomial
    -- * References
    -- $references
  ) where

import Control.Arrow           (second,(***))

import Data.Word
import Data.Int
import Data.Bits
import qualified Data.Vector.Generic         as G
import           Data.Vector.Generic           ((++))
import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Unboxed         as U
import qualified Data.Vector                 as V
import Data.Vector.Generic (Vector)
import Numeric.SpecFunctions (logFactorial)
import System.Random.Stateful

import Prelude hiding ((++))

-- | A lookup table for arbitrary discrete distributions. It allows
-- the generation of random variates in /O(1)/. Note that probability
-- is quantized in units of @1/2^32@, and all distributions with
-- infinite support (e.g. Poisson) should be truncated.
data CondensedTable v a =
  CondensedTable
  {-# UNPACK #-} !Word64 !(v a) -- Lookup limit and first table
  {-# UNPACK #-} !Word64 !(v a) -- Second table
  {-# UNPACK #-} !Word64 !(v a) -- Third table
  !(v a)                        -- Last table

-- Implementation note. We have to store lookup limit in Word64 since
-- we need to accomodate two cases. First is when we have no values in
-- lookup table, second is when all elements are there
--
-- Both are pretty easy to realize. For first one probability of every
-- outcome should be less then 1/256, latter arise when probabilities
-- of two outcomes are [0.5,0.5]

-- | A 'CondensedTable' that uses unboxed vectors.
type CondensedTableU = CondensedTable U.Vector

-- | A 'CondensedTable' that uses boxed vectors, and is able to hold
-- any type of element.
type CondensedTableV = CondensedTable V.Vector

-- | Generate a random value using a condensed table.
genFromTable :: (StatefulGen g m, Vector v a) => CondensedTable v a -> g -> m a
{-# INLINE genFromTable #-}
genFromTable table gen = do
  w <- uniformM gen
  return $ lookupTable table $ fromIntegral (w :: Word32)

lookupTable :: Vector v a => CondensedTable v a -> Word64 -> a
{-# INLINE lookupTable #-}
lookupTable (CondensedTable na aa nb bb nc cc dd) i
  | i < na    = aa `at` ( i       `shiftR` 24)
  | i < nb    = bb `at` ((i - na) `shiftR` 16)
  | i < nc    = cc `at` ((i - nb) `shiftR` 8 )
  | otherwise = dd `at` ( i - nc)
  where
    at arr j = G.unsafeIndex arr (fromIntegral j)

----------------------------------------------------------------
-- Table generation
----------------------------------------------------------------

-- | Generate a condensed lookup table from a list of outcomes with
-- given probabilities. The vector should be non-empty and the
-- probabilites should be non-negative and sum to 1. If this is not
-- the case, this algorithm will construct a table for some
-- distribution that may bear no resemblance to what you intended.
tableFromProbabilities
    :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32)
       => v (a, Double) -> CondensedTable v a
{-# INLINE tableFromProbabilities #-}
tableFromProbabilities v
  | G.null tbl = pkgError "tableFromProbabilities" "empty vector of outcomes"
  | otherwise  = tableFromIntWeights $ G.map (second $ toWeight . (* mlt)) tbl
  where
    -- 2^32. N.B. This number is exatly representable.
    mlt = 4.294967296e9
    -- Drop non-positive probabilities
    tbl = G.filter ((> 0) . snd) v
    -- Convert Double weight to Word32 and avoid overflow at the same
    -- time. It's especially dangerous if one probability is
    -- approximately 1 and others are 0.
    toWeight w | w > mlt - 1 = 2^(32::Int) - 1
               | otherwise   = round w

-- | Same as 'tableFromProbabilities' but treats number as weights not
-- probilities. Non-positive weights are discarded, and those
-- remaining are normalized to 1.
tableFromWeights
    :: (Vector v (a,Word32), Vector v (a,Double), Vector v a, Vector v Word32)
       => v (a, Double) -> CondensedTable v a
{-# INLINE tableFromWeights #-}
tableFromWeights = tableFromProbabilities . normalize . G.filter ((> 0) . snd)
  where
    normalize v
      | G.null v  = pkgError "tableFromWeights" "no positive weights"
      | otherwise = G.map (second (/ s)) v
      where
        -- Explicit fold is to avoid 'Vector v Double' constraint
        s = G.foldl' (flip $ (+) . snd) 0 v

-- | Generate a condensed lookup table from integer weights. Weights
-- should sum to @2^32@ at least approximately. This function will
-- correct small deviations from @2^32@ such as arising from rounding
-- errors. But for large deviations it's likely to product incorrect
-- result with terrible performance.
tableFromIntWeights :: (Vector v (a,Word32), Vector v a, Vector v Word32)
                    => v (a, Word32)
                    -> CondensedTable v a
{-# INLINE tableFromIntWeights #-}
tableFromIntWeights v
  | n == 0    = pkgError "tableFromIntWeights" "empty table"
    -- Single element tables should be treated sepately. Otherwise
    -- they will confuse correctWeights
  | n == 1    = let m = 2^(32::Int) - 1 -- Works for both Word32 & Word64
                in CondensedTable
                   m (G.replicate 256 $ fst $ G.head tbl)
                   m  G.empty
                   m  G.empty
                      G.empty
  | otherwise = CondensedTable
                na aa
                nb bb
                nc cc
                   dd
  where
    -- We must filter out zero-probability outcomes because they may
    -- confuse weight correction algorithm
    tbl   = G.filter ((/=0) . snd) v
    n     = G.length tbl
    -- Corrected table
    table = uncurry G.zip $ id *** correctWeights $ G.unzip tbl
    -- Make condensed table
    mkTable  d =
      G.concatMap (\(x,w) -> G.replicate (fromIntegral $ digit d w) x) table
    len = fromIntegral . G.length
    -- Tables
    aa = mkTable 0
    bb = mkTable 1
    cc = mkTable 2
    dd = mkTable 3
    -- Offsets
    na =       len aa `shiftL` 24
    nb = na + (len bb `shiftL` 16)
    nc = nb + (len cc `shiftL` 8)

-- Calculate N'th digit base 256
digit :: Int -> Word32 -> Word32
digit 0 x =  x `shiftR` 24
digit 1 x = (x `shiftR` 16) .&. 0xff
digit 2 x = (x `shiftR` 8 ) .&. 0xff
digit 3 x =  x .&. 0xff
digit _ _ = pkgError "digit" "the impossible happened!?"
{-# INLINE digit #-}

-- Correct integer weights so they sum up to 2^32. Array of weight
-- should contain at least 2 elements.
correctWeights :: G.Vector v Word32 => v Word32 -> v Word32
{-# INLINE correctWeights #-}
correctWeights v = G.create $ do
  let
    -- Sum of weights
    s = G.foldl' (flip $ (+) . fromIntegral) 0 v :: Int64
    -- Array size
    n = G.length v
  arr <- G.thaw v
  -- On first pass over array adjust only entries which are larger
  -- than `lim'. On second and subsequent passes `lim' is set to 1.
  --
  -- It's possibly to make this algorithm loop endlessly if all
  -- weights are 1 or 0.
  let loop lim i delta
        | delta == 0 = return ()
        | i >= n     = loop 1 0 delta
        | otherwise  = do
            w <- M.read arr i
            case () of
              _| w < lim   -> loop lim (i+1) delta
               | delta < 0 -> M.write arr i (w + 1) >> loop lim (i+1) (delta + 1)
               | otherwise -> M.write arr i (w - 1) >> loop lim (i+1) (delta - 1)
  loop 255 0 (s - 2^(32::Int))
  return arr

-- | Create a lookup table for the Poisson distibution. Note that
-- table construction may have significant cost. For λ < 100 it
-- takes as much time to build table as generation of 1000-30000
-- variates.
tablePoisson :: Double -> CondensedTableU Int
tablePoisson = tableFromProbabilities . make
  where
    make lam
      | lam < 0    = pkgError "tablePoisson" "negative lambda"
      | lam < 22.8 = U.unfoldr unfoldForward (exp (-lam), 0)
      | otherwise  = U.unfoldr unfoldForward (pMax, nMax)
                  ++ U.tail (U.unfoldr unfoldBackward (pMax, nMax))
      where
        -- Number with highest probability and its probability
        --
        -- FIXME: this is not ideal precision-wise. Check if code
        --        from statistics gives better precision.
        nMax = floor lam :: Int
        pMax = exp $ fromIntegral nMax * log lam - lam - logFactorial nMax
        -- Build probability list
        unfoldForward (p,i)
          | p < minP  = Nothing
          | otherwise = Just ( (i,p)
                             , (p * lam / fromIntegral (i+1), i+1)
                             )
        -- Go down
        unfoldBackward (p,i)
          | p < minP  = Nothing
          | otherwise = Just ( (i,p)
                             , (p / lam * fromIntegral i, i-1)
                             )
    -- Minimal representable probability for condensed tables
    minP = 1.1641532182693481e-10 -- 2**(-33)

-- | Create a lookup table for the binomial distribution.
tableBinomial :: Int            -- ^ Number of tries
              -> Double         -- ^ Probability of success
              -> CondensedTableU Int
tableBinomial n p = tableFromProbabilities makeBinom
  where
  makeBinom
    | n <= 0         = pkgError "tableBinomial" "non-positive number of tries"
    | p == 0         = U.singleton (0,1)
    | p == 1         = U.singleton (n,1)
    | p > 0 && p < 1 = U.unfoldrN (n + 1) unfolder ((1-p)^n, 0)
    | otherwise      = pkgError "tableBinomial" "probability is out of range"
    where
      h = p / (1 - p)
      unfolder (t,i) = Just ( (i,t)
                            , (t * (fromIntegral $ n + 1 - i1) * h / fromIntegral i1, i1) )
        where i1 = i + 1

pkgError :: String -> String -> a
pkgError func err =
    error . concat $ ["System.Random.MWC.CondensedTable.", func, ": ", err]

-- $references
--
-- * Wang, J.; Tsang, W. W.; G. Marsaglia (2004), Fast Generation of
--   Discrete Random Variables, /Journal of Statistical Software,
--   American Statistical Association/, vol. 11(i03).
--   <http://ideas.repec.org/a/jss/jstsof/11i03.html>

max_stars_repo_name: Shimuuar/mwc-random

max_issues_repo_name: Shimuuar/mwc-random

max_forks_repo_name: Shimuuar/mwc-random


text:
module GaussianPinwheel where

import           Control.Parallel.Strategies
import           Data.Array.Repa                          as R
import           Data.Complex                             as C
import           Data.List                                as L
import           Data.Vector.Storable                     as VS
import           Data.Vector.Unboxed                      as VU
import           STC.Plan
import           Filter.Utils
import qualified FourierPinwheel.Filtering                as FP
import           FourierPinwheel.GaussianEnvelopePinwheel
import           Graphics.Rendering.Chart.Backend.Cairo
import           Graphics.Rendering.Chart.Easy
import           Image.IO
import           Math.Gamma
import           Pinwheel.FourierSeries2D
import           System.Directory
import           System.Environment
import           System.FilePath
import           Text.Printf
import           Utils.BLAS
import           Utils.Distribution
import           Utils.List

main = do
  args@(numR2FreqsStr:numThetaFreqStr:numRFreqStr:numPointStr:numThetaStr:numRStr:angularFreqStr:radialFreqStr:stdStr:sigmaStr:_) <-
    getArgs
  let numR2Freqs = read numR2FreqsStr :: Int
      numThetaFreq = read numThetaFreqStr :: Int
      numRFreq = read numRFreqStr :: Int
      numPoint = read numPointStr :: Int
      numTheta = read numThetaStr :: Int
      numR = read numRStr :: Int
      angularFreq = read angularFreqStr :: Int
      radialFreq = read radialFreqStr :: Int
      std = read stdStr :: Double
      sigma = read sigmaStr :: Double
      periodR2 = fromIntegral numR2Freqs
      periodEnv = periodR2 ^ 2 / 4
      folderPath = "output/test/GaussianPinwheel"
  createDirectoryIfMissing True folderPath
  --DFT Plan
  initVec <- generateRadomVector (numR2Freqs ^ 2)
  lock <- getFFTWLock
  plan <-
    fst <$>
    (idft1dGPlan lock emptyPlan [numR2Freqs, numR2Freqs] [0, 1] initVec >>= \(plan, vec) ->
       dft1dGPlan lock plan [numR2Freqs, numR2Freqs] [0, 1] vec)
  let a = 1 / (std * sqrt 2)
      pinwheelCoef =
        centerHollowArray numR2Freqs $
        createFrequencyArray
          numR2Freqs
          (gaussianPinwheelFourierCoefficients
             numR2Freqs
             periodR2
             a
             sigma
             angularFreq
             radialFreq
             periodEnv)
  plotImageRepaComplex (folderPath </> "GaussianPinwheelCoef.png") .
    ImageRepa 8 . computeS . extend (Z :. (1 :: Int) :. All :. All) $
    pinwheelCoef
  pinwheel <-
    dftExecute plan (DFTPlanID IDFT1DG [numR2Freqs, numR2Freqs] [0, 1]) .
    VU.convert . toUnboxed . computeS . makeFilter2D $
    pinwheelCoef
  plotImageRepaComplex (folderPath </> "GaussianPinwheel.png") .
    ImageRepa 8 .
    computeS .
    extend (Z :. (1 :: Int) :. All :. All) .
    makeFilter2DInverse .
    fromUnboxed (Z :. numR2Freqs :. numR2Freqs) . VS.convert $
    pinwheel
  let pinwheel1 =
        computeS $
        createFrequencyArray
          numR2Freqs
          (\phi rho ->
             ((exp ((-1) * (rho / std) ^ 2 / 2)) :+ 0) *
             (rho :+ 0) **
             (sigma :+ ((-2) * pi * fromIntegral radialFreq / log periodEnv)) *
             cis ((-phi) * fromIntegral angularFreq))
  pinwheelCoef1 <-
    dftExecute plan (DFTPlanID DFT1DG [numR2Freqs, numR2Freqs] [0, 1]) .
    VU.convert . toUnboxed . computeS . makeFilter2D $
    pinwheel1
  plotImageRepaComplex (folderPath </> "GaussianPinwheelCoef1.png") .
    ImageRepa 8 .
    computeS .
    extend (Z :. (1 :: Int) :. All :. All) .
    makeFilter2DInverse .
    fromUnboxed (Z :. numR2Freqs :. numR2Freqs) . VS.convert $
    pinwheelCoef1
  plotImageRepaComplex (folderPath </> "GaussianPinwheel1.png") .
    ImageRepa 8 .
    computeS .
    extend (Z :. (1 :: Int) :. All :. All) .
    -- makeFilter2DInverse .
    fromUnboxed (Z :. numR2Freqs :. numR2Freqs) . VS.convert . toUnboxed $
    pinwheel1

max_stars_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_issues_repo_name: XinhuaZhang/Stochastic-Completion-Field

max_forks_repo_name: XinhuaZhang/Stochastic-Completion-Field


text:
import Data.Complex (cis, phase)

meanAngle = (/ pi) . (* 180) . phase . sum . map (cis . (/ 180) . (* pi))

main = mapM_ (\angles -> putStrLn $ "The mean angle of " ++ show angles ++ " is: " ++ show (meanAngle angles) ++ " degrees")
       [[350, 10], [90, 180, 270, 360], [10, 20, 30]]

max_stars_repo_name: djgoku/RosettaCodeData

max_issues_repo_name: djgoku/RosettaCodeData

max_forks_repo_name: djgoku/RosettaCodeData