boltzmann-samplers (empty) → 0.1.0.0

raw patch · 14 files changed

+1923/−0 lines, 14 filesdep +MonadRandomdep +QuickCheckdep +adsetup-changed

Dependencies added: MonadRandom, QuickCheck, ad, base, boltzmann-samplers, containers, criterion, deepseq, hashable, hmatrix, ieee754, mtl, optparse-generic, testing-feat, transformers, unordered-containers, vector

Files

• LICENSE +20/−0
• README.md +54/−0
• Setup.hs +2/−0
• bench/binaryTree.hs +100/−0
• boltzmann-samplers.cabal +90/−0
• src/Boltzmann/Data.hs +313/−0
• src/Boltzmann/Data/Common.hs +39/−0
• src/Boltzmann/Data/Data.hs +152/−0
• src/Boltzmann/Data/Oracle.hs +506/−0
• src/Boltzmann/Data/Types.hs +197/−0
• src/Boltzmann/Solver.hs +69/−0
• src/Boltzmann/Species.hs +220/−0
• test/Test/Stats.hs +76/−0
• test/tree.hs +85/−0

@@ -0,0 +1,20 @@+Copyright 2017 Li-yao Xia+Copyright 2017 Li-yao Xia++Permission is hereby granted, free of charge, to any person obtaining a copy of+this software and associated documentation files (the "Software"), to deal in+the Software without restriction, including without limitation the rights to+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies+of the Software, and to permit persons to whom the Software is furnished to do+so, subject to the following conditions:++The above copyright notice and this permission notice shall be included in all+copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE+SOFTWARE.

@@ -0,0 +1,54 @@+Boltzmann samplers [![Hackage](https://img.shields.io/hackage/v/boltzmann-samplers.svg)](https://hackage.haskell.org/package/generic-random) [![Build Status](https://travis-ci.org/Lysxia/generic-random.svg)](https://travis-ci.org/Lysxia/boltzmann-samplers)+==================++`Boltzmann.Data`+----------------++Define sized random generators for `Data.Data` generic types.++```haskell+    {-# LANGUAGE DeriveDataTypeable #-}++    import Data.Data+    import Test.QuickCheck+    import Boltzmann.Data++    data Term = Lambda Int Term | App Term Term | Var Int+      deriving (Show, Data)++    instance Arbitrary Term where+      arbitrary = sized $ generatorPWith [positiveInts]++    positiveInts :: Alias Gen+    positiveInts =+      alias $ \() -> fmap getPositive arbitrary :: Gen Int++    main = sample (arbitrary :: Gen Term)+```++- Objects of the same size (number of constructors) occur with the same+  probability (see Duchon et al., references below).+- Implements rejection sampling and pointing.+- Works with QuickCheck and MonadRandom, but also similar user-defined monads+  for randomness (just implement `MonadRandomLike`).+- Can be tweaked somewhat with user defined generators.++`Boltzmann.Species`+-------------------++An experimental interface to obtain Boltzmann samplers from an applicative+specification of a combinatorial system.++No documentation (yet).++References+----------++- The core theory of Boltzmann samplers is described in+  [Boltzmann Samplers for the Random Generation of Combinatorial Structures](http://algo.inria.fr/flajolet/Publications/DuFlLoSc04.pdf),+  P. Duchon, P. Flajolet, G. Louchard, G. Schaeffer.++- The numerical evaluation of recursively defined generating functions+  is taken from+  [Boltzmann Oracle for Combinatorial Systems](http://www.dmtcs.org/pdfpapers/dmAI0132.pdf),+  C. Pivoteau, B. Salvy, M. Soria.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,100 @@+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE TemplateHaskell #-}++module Main where++import Control.Applicative+import Control.Monad+import Control.Monad.Trans.Class+import Data.Bool+import Data.Data+import Data.Functor+import GHC.Generics+import Control.DeepSeq+import Criterion.Main+import Test.Feat+import Test.QuickCheck+import Test.QuickCheck.Gen+import Test.QuickCheck.Random+import Control.Exception ( evaluate )++import Boltzmann.Data+import Boltzmann.Data.Data+import Boltzmann.Data.Types++data T = N T T | L+  deriving (Eq, Ord, Show, Data, Typeable, Generic)++instance NFData T++deriveEnumerable ''T++size :: Num a => T -> a+size L = 1+size (N l r) = 1 + size l + size r++gen1 :: Int -> Gen T+gen1 n = runRejectT (tolerance epsilon (n + 1)) gen'+  where+    gen' = incr >> lift arbitrary >>= bool (return L) (liftA2 N gen' gen')++gen2 :: Int -> Gen T+gen2 n = g+  where+    (minSize, maxSize) = tolerance epsilon (n + 1)+    g = gen' 0 (\m t -> if m < minSize then g else return t)+    gen' n k | n >= maxSize = g+    gen' n k =+      arbitrary >>= bool+        (k (n+1) L)+        (gen' (n+1) $ \m l -> gen' m $ \m r -> k m (N l r))++genFeat :: Int -> Gen T+genFeat = uniform++main = newQCGen >>= \g -> defaultMain $ liftA2 (\n f -> f n g)+  [4 ^ e | e <- [1 .. 6]]++  -- Singular rejection sampling+  [ bg "handwritten1" gen1+  , bg "handwritten2" gen2++  , bg "feat" genFeat++  -- Pointed generator+  , bg "P" generatorP'++  -- Pointed generator with rejection sampling+  , bg "PR" generatorPR'++  , bg "SR" generatorSR++  -- Sized rejection sampling+  , bg "R" generatorR'++  -- Sized rejection sampling, not memoizing oracle+  , bg' "R-recomp" generatorR'++  -- Pointed generator, not memoizing oracle+  , bg' "P-recomp" generatorP'+  ]++bg, bg' :: String -> (Int -> Gen T) -> Int -> QCGen -> Benchmark+bg name gen n g =+  bench (name ++ "_" ++ show n) $ nf f g+  where+    go 0 = return (0 :: Int)+    go k = liftA2 (\t s -> size t + s) gg (go (k-1))+    gg = gen n+    f g = unGen (go 100) g 0++bg' name gen n g =+  bench (name ++ "_" ++ show n) $ nf f (n, g)+  where+    go n 0 = return (0 :: Int)+    go n k = liftA2 (\t s -> size t + s) (gen n) (go n (k-1))+    f (n, g) = unGen (go n 100) g 0++avgSize :: [T] -> Double+avgSize ts = sum (fmap size ts) / fromIntegral (length ts)

@@ -0,0 +1,90 @@+name:                boltzmann-samplers+version:             0.1.0.0+synopsis:            Uniform random generators+description:++  Random generators with a uniform distribution conditioned+  to a given size.++  See also @<http://hackage.haskell.org/package/testing-feat testing-feat>@,+  which is currently a faster method with similar results.++homepage:            https://github.com/Lysxia/boltzmann-samplers#readme+license:             MIT+license-file:        LICENSE+author:              Li-yao Xia+maintainer:          lysxia@gmail.com+category:            Data, Generic, Random+build-type:          Simple+extra-source-files:  README.md+cabal-version:       >=1.10++flag test+  Description:+    Enable testing. Disabled by default because the current test suite+    is slow and can fail with non-zero probability.+  Manual:  True+  Default: False++library+  hs-source-dirs:      src+  exposed-modules:+    Boltzmann.Data+    Boltzmann.Data.Data+    Boltzmann.Data.Common+    Boltzmann.Data.Oracle+    Boltzmann.Data.Types+    Boltzmann.Solver+    Boltzmann.Species+  build-depends:+    ad,+    base >= 4.9 && < 5,+    containers,+    hashable,+    hmatrix,+    ieee754,+    unordered-containers,+    MonadRandom,+    mtl,+    QuickCheck,+    transformers,+    vector+  default-language:    Haskell2010++source-repository head+  type:     git+  location: https://github.com/Lysxia/boltzmann-samplers++test-suite test-tree+  type:             exitcode-stdio-1.0+  hs-source-dirs:   test+  main-is:          tree.hs+  default-language: Haskell2010+  other-modules:+    Test.Stats+  if flag(test)+    build-depends:+      base,+      QuickCheck,+      optparse-generic,+      boltzmann-samplers+  else+    buildable: False++benchmark bench-binarytree+  type:             exitcode-stdio-1.0+  hs-source-dirs:   bench+  main-is:          binaryTree.hs+  default-language: Haskell2010+  ghc-options: -O2+  if flag(test)+    build-depends:+      base,+      criterion,+      deepseq,+      QuickCheck,+      transformers,+      testing-feat,+      boltzmann-samplers+  else+    buildable: False

@@ -0,0 +1,313 @@+-- | Generic Boltzmann samplers.+--+-- Here, the words "/sampler/" and "/generator/" are used interchangeably.+--+-- Given an algebraic datatype:+--+-- > data A = A1 B C | A2 D+--+-- a Boltzmann sampler is recursively defined by choosing a constructor with+-- some fixed distribution, and /independently/ generating values for the+-- corresponding fields with the same method.+--+-- A key component is the aforementioned distribution, defined for every type+-- such that the resulting generator produces a finite value in the end. These+-- distributions are obtained from a precomputed object called /oracle/, which+-- we will not describe further here.+--+-- Oracles depend on the target size of the generated data (except for singular+-- samplers), and can be fairly expensive to compute repeatedly, hence some of+-- the functions below attempt to avoid (re)computing too many of them even+-- when the required size changes.+--+-- When these functions are specialized, oracles are memoized and will be+-- reused for different sizes.++module Boltzmann.Data (+  Size',+  -- * Main functions+  -- $sized+  generatorSR,+  generatorP,+  generatorPR,+  generatorR,+  -- ** Fixed size+  -- $fixed+  generatorP',+  generatorPR',+  generatorR',+  generator',+  -- * Generators with aliases+  -- $aliases+  generatorSRWith,+  generatorPWith,+  generatorPRWith,+  generatorRWith,+  -- ** Fixed size+  generatorPWith',+  generatorPRWith',+  generatorRWith',+  generatorWith',+  -- * Other generators+  -- $other+  Points,+  generatorM,+  generatorMR,+  generator_,+  generatorR_,+  -- * Auxiliary definitions+  -- ** Type classes+  MonadRandomLike (..),+  AMonadRandom (..),+  -- ** Alias+  alias,+  aliasR,+  coerceAlias,+  coerceAliases,+  Alias (..),+  AliasR,+  ) where++import Data.Data+import Boltzmann.Data.Data+import Boltzmann.Data.Types++-- * Main functions++-- $sized+--+-- === Suffixes+--+-- [@S@] Singular sampler.+--+--     This works with recursive tree-like structures, as opposed to (lists of)+--     structures with bounded size. More precisely, the generating function of+--     the given type should have a finite radius of convergence, with a+--     singularity of a certain kind (see Duchon et al., reference in the+--     README), so that the oracle can be evaluated at that point.+--+--     This has the advantage of using the same oracle for all size parameters,+--     which simply specify a target size interval.+--+-- [@P@] Generator of pointed values.+--+--     It usually has a flatter distribution of sizes than a simple Boltzmann+--     sampler, making it an efficient alternative to rejection sampling.+--+--     It also works on more types, particularly lists and finite types,+--     but relies on multiple oracles.+--+-- [@R@] Rejection sampling.+--+--     These generators filter out values whose sizes are not within some+--     interval. In the first two sections, that interval is implicit:+--     @[(1-'epsilon')*size', (1+'epsilon')*size']@, for @'epsilon' = 0.1@.+--+--     The generator restarts as soon as it has produced more constructors than+--     the upper bound, this strategy is called /ceiled rejection sampling/.+--+-- = Pointing+--+-- The /pointing/ of a type @t@ is a derived type whose values are essentially+-- values of type @t@, with one of their constructors being "pointed".+-- Alternatively, we may turn every constructor into variants that indicate+-- the position of points.+--+-- @+--   -- Original type+--   data Tree = Node Tree Tree | Leaf+--   -- Pointing of Tree+--   data Tree'+--     = Tree' Tree -- Point at the root+--     | Node'0 Tree' Tree -- Point to the left+--     | Node'1 Tree Tree' -- Point to the right+-- @+--+-- Pointed values are easily mapped back to the original type by erasing the+-- point. Pointing makes larger values occur much more frequently, while+-- preserving the uniformness of the distribution conditionally to a fixed+-- size.+--++-- | @+--   'generatorSR' :: Int -> 'Gen' a+--   'asMonadRandom' . 'generatorSR' :: 'MonadRandom' m => Int -> m a+-- @+--+-- Singular ceiled rejection sampler.+generatorSR :: (Data a, MonadRandomLike m) => Size' -> m a+generatorSR = generatorSRWith []++-- | @+--   'generatorP' :: Int -> 'Gen' a+--   'asMonadRandom' . 'generatorP' :: 'MonadRandom' m => Int -> m a+-- @+--+-- Generator of pointed values.++generatorP :: (Data a, MonadRandomLike m) => Size' -> m a+generatorP = generatorPWith []++-- | Pointed generator with rejection.+generatorPR :: (Data a, MonadRandomLike m) => Size' -> m a+generatorPR = generatorPRWith []++-- | Generator with rejection and dynamic average size.+generatorR :: (Data a, MonadRandomLike m) => Size' -> m a+generatorR = generatorRWith []++-- ** Fixed size++-- $fixed+-- The @'@ suffix indicates functions which do not do any+-- precomputation before passing the size parameter.+--+-- This means that oracles are computed from scratch for every size value,+-- which may incur a significant overhead.++-- | Pointed generator.+generatorP' :: (Data a, MonadRandomLike m) => Size' -> m a+generatorP' = generatorPWith' []++-- | Pointed generator with rejection.+generatorPR' :: (Data a, MonadRandomLike m) => Size' -> m a+generatorPR' = generatorPRWith' []++-- | Ceiled rejection sampler with given average size.+generatorR' :: (Data a, MonadRandomLike m) => Size' -> m a+generatorR' = generatorRWith' []++-- | Basic boltzmann sampler with no optimization.+generator' :: (Data a, MonadRandomLike m) => Size' -> m a+generator' = generatorWith' []++-- * Generators with aliases++-- $aliases+-- Boltzmann samplers can normally be defined only for types @a@ such that:+--+-- - they are instances of 'Data';+-- - the set of types of subterms of values of type @a@ is finite;+-- - and all of these types have at least one finite value (i.e., values with+--   finitely many constructors).+--+-- Examples of misbehaving types are:+--+-- - @a -> b -- Not Data@+-- - @data E a = L a | R (E [a]) -- Contains a, [a], [[a]], [[[a]]], etc.@+-- - @data I = C I -- No finite value@+--+-- = Alias+--+-- The 'Alias' type works around these limitations ('AliasR' for rejection+-- samplers).+-- This existential wrapper around a user-defined function @f :: a -> m b@+-- makes @boltzmann-samplers@ view occurences of the type @b@ as @a@ when+-- processing a recursive system of types, possibly stopping some infinite+-- unrolling of type definitions. When a value of type @b@ needs to be+-- generated, it generates an @a@ which is passed to @f@.+--+-- @+--   let+--     as = ['aliasR' $ \\() -> return (L []) :: 'Gen' (E [[Int]])]+--   in+--     'generatorSRWith' as 'asGen' :: 'Size' -> 'Gen' (E Int)+-- @+--+-- Another use case is to plug in user-defined generators where the default is+-- not satisfactory, for example, to generate positive @Int@s:+--+-- @+--   let+--     as = ['alias' $ \\() -> 'choose' (0, 100) :: 'Gen' Int)]+--   in+--     'generatorPWith' as 'asGen' :: 'Size' -> 'Gen' [Int]+-- @+--+-- or to modify the weights assigned to some types. In particular, in some+-- cases it seems preferable to make @String@ (and @Text@) have the same weight+-- as @Int@ and @()@.+--+-- @+--   let+--     as = ['alias' $ \\() -> arbitrary :: 'Gen' String]+--   in+--     'generatorPWith' as 'asGen' :: 'Size' -> 'Gen' (Either Int String)+-- @++generatorSRWith+  :: (Data a, MonadRandomLike m) => [AliasR m] -> Size' -> m a+generatorSRWith aliases =+  generatorR_ aliases 0 Nothing . tolerance epsilon++generatorPRWith+  :: (Data a, MonadRandomLike m) => [AliasR m] -> Size' -> m a+generatorPRWith aliases size' =+  generatorMR aliases 1 size' (tolerance epsilon size')++generatorPWith+  :: (Data a, MonadRandomLike m) => [Alias m] -> Size' -> m a+generatorPWith aliases = generatorM aliases 1++generatorRWith+  :: (Data a, MonadRandomLike m) => [AliasR m] -> Size' -> m a+generatorRWith aliases size' =+  generatorMR aliases 0 size' (tolerance epsilon size')++-- ** Fixed size++generatorPWith'+  :: (Data a, MonadRandomLike m) => [Alias m] -> Size' -> m a+generatorPWith' aliases = generator_ aliases 1 . Just++generatorPRWith'+  :: (Data a, MonadRandomLike m) => [AliasR m] -> Size' -> m a+generatorPRWith' aliases size' =+  generatorR_ aliases 1 (Just size') (tolerance epsilon size')++generatorRWith'+  :: (Data a, MonadRandomLike m) => [AliasR m] -> Size' -> m a+generatorRWith' aliases size' =+  generatorR_ aliases 0 (Just size') (tolerance epsilon size')++generatorWith'+  :: (Data a, MonadRandomLike m) => [Alias m] -> Size' -> m a+generatorWith' aliases = generator_ aliases 0 . Just++-- * Other generators++-- $other Used in the implementation of the generators above.+-- These also allow to apply pointing more than once.+--+-- === Suffixes+--+-- [@M@] Sized generators are memoized for some sparsely chosen values of+-- sizes. Subsequently supplied sizes are approximated by the closest larger+-- value. This strategy avoids recomputing too many oracles. Aside from+-- singular samplers, all other generators above not marked by @'@ use this.+--+-- [@_@] If the size parameter is @Nothing@, produces the singular generator+-- (associated with the suffix @S@); otherwise the generator produces values+-- with average size equal to the given value.++generatorM+  :: (Data a, MonadRandomLike m)+  => [Alias m] -> Points -> Size' -> m a+generatorM = memo make apply++generatorMR+  :: (Data a, MonadRandomLike m)+  => [AliasR m] -> Points -> Size' -> (Size', Size') -> m a+generatorMR = memo makeR applyR++-- | Boltzmann sampler without rejection.+generator_+  :: (Data a, MonadRandomLike m)+  => [Alias m] -> Points -> Maybe Size' -> m a+generator_ aliases = apply (make aliases [])++-- | Boltzmann sampler with rejection.+generatorR_+  :: (Data a, MonadRandomLike m)+  => [AliasR m] -> Points -> Maybe Size' -> (Size', Size') -> m a+generatorR_ aliases = applyR (makeR aliases [])

@@ -0,0 +1,39 @@+-- | General helper functions++module Boltzmann.Data.Common where++frequencyWith+  :: (Ord r, Num r, Monad m) => (r -> m r) -> [(r, m a)] -> m a+frequencyWith _ [(_, a)] = a+frequencyWith randomR as = randomR total >>= select as+  where+    total = (sum . fmap fst) as+    select ((w, a) : as) x+      | x < w = a+      | otherwise = select as (x - w)+    select _ _ = (snd . head) as+    -- That should not happen in theory, but floating point might be funny.++-- | @partitions k n@: lists of non-negative integers of length @n@ with sum+-- less than or equal to @k@.+partitions :: Int -> Int -> [[Int]]+partitions _ 0 = [[]]+partitions k n = do+  p <- [0 .. k]+  (p :) <$> partitions (k - p) (n - 1)++-- | Binomial coefficient.+--+-- > binomial n k == factorial n `div` (factorial k * factorial (n-k))+binomial :: Int -> Int -> Integer+binomial = \n k -> pascal !! n !! k+  where+    pascal = [1] : fmap nextRow pascal+    nextRow r = zipWith (+) (0 : r) (r ++ [0])++-- | Multinomial coefficient.+--+-- > multinomial n ps == factorial n `div` product [factorial p | p <- ps]+multinomial :: Int -> [Int] -> Integer+multinomial _ [] = 1+multinomial n (p : ps) = binomial n p * multinomial (n - p) ps

@@ -0,0 +1,152 @@+-- | Internal module+--+-- Derive Boltzmann samplers for SYB-generic types.++{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE RecordWildCards #-}++module Boltzmann.Data.Data where++import Control.Arrow ( (&&&) )+import Control.Applicative+import Data.Data+import Data.Foldable+import Data.Maybe+import qualified Data.HashMap.Lazy as HashMap+import Boltzmann.Data.Oracle+import Boltzmann.Data.Types++-- | Sized generator.+data SG r = SG+  { minSize :: Size+  , maxSizeM :: Maybe Size+  , runSG :: Points -> Maybe Double -> r+  , runSmallG :: Points -> r+  } deriving Functor++-- | Number of pointing iterations.+type Points = Int++rangeSG :: SG r -> (Size, Maybe Size)+rangeSG = minSize &&& maxSizeM++-- | For documentation.+applySG :: SG r -> Points -> Maybe Double -> r+applySG SG{..} k sizeM+  | Just minSize == maxSizeM = runSG k (fmap fromIntegral maxSizeM)+  | Just size <- sizeM, size <= fromIntegral minSize =+      error "Target size too small."+  | Just True <- liftA2 ((<=) . fromIntegral) maxSizeM sizeM =+      error "Target size too large."+  | Nothing <- sizeM, Just _ <- maxSizeM =+      error "Cannot make singular sampler for finite type."+  | otherwise = runSG k sizeM++-- * Helper functions++make :: (Data a, MonadRandomLike m)+  => [Alias m] -> proxy a -> SG (m a)+make aliases a =+  SG minSize maxSizeM make' makeSmall+  where+    dd = collectTypes aliases a+    t = typeRep a+    i = case index dd #! t of+      Left j -> fst (xedni' dd #! j)+      Right i -> i+    minSize = natToInt $ fst (lTerm dd #! i)+    maxSizeM = HashMap.lookup i (degree dd)+    make' k sizeM = getGenerator dd' generators a k+      where+        dd' = dds !! k+        oracle = makeOracle dd' t sizeM+        generators = makeGenerators dd' oracle+    makeSmall k = getSmallGenerator dd' (smallGenerators dd') a+      where+        dd' = dds !! k+    dds = iterate point dd++makeR :: (Data a, MonadRandomLike m)+  => [AliasR m] -> proxy a+  -> SG ((Size, Size) -> m a)+makeR aliases a = fmap (flip runRejectT) (make aliases a)++-- | The size of a value is its number of constructors.+--+-- Here, however, the 'Size'' type is interpreted differently to make better+-- use of QuickCheck's size parameter provided by the 'Test.QuickCheck.sized'+-- combinator, so that we generate non-trivial data even at very small size+-- values.+--+-- For infinite types, with objects of unbounded sizes @> minSize@, given a+-- parameter @delta :: 'Size''@, the produced values have an average size close+-- to @minSize + delta@.+--+-- For example, values of type @Either () [Bool]@ have at least two constructors,+-- so+--+-- @+--   'generator' delta :: 'Gen' (Either () [Bool])+-- @+--+-- will target sizes close to @2 + delta@;+-- the offset becomes less noticeable as @delta@ grows to infinity.+--+-- For finite types with sizes in @[minSize, maxSize]@, the target expected+-- size is obtained by clamping a 'Size'' to @[0, 99]@ and applying an affine+-- mapping.+type Size' = Int++rescale :: SG r -> Size' -> Double+rescale (SG minSize (Just maxSize) _ _) size' =+  fromIntegral minSize + fromIntegral (min 99 size' * (maxSize - minSize)) / 100+rescale (SG minSize Nothing _ _) size' = fromIntegral (minSize + size')++apply :: SG r -> Points -> Maybe Size' -> r+apply sg k (Just 0) = runSmallG sg k+apply sg k size' = runSG sg k (fmap (rescale sg) size')++applyR :: SG ((Size, Size) -> r) -> Points -> Maybe Size' -> (Size', Size') -> r+applyR sg k size' = apply sg k size' . rescaleInterval sg++rescaleInterval :: SG r -> (Size', Size') -> (Size, Size)+rescaleInterval sg (a', b') = (a, b)+  where+    a = (clamp . floor .rescale sg) a'+    b = (clamp . ceiling . rescale sg) b'+    clamp x+      | Just maxSize <- maxSizeM sg, x >= 100 = maxSize+      | otherwise = x++-- | > 'epsilon' = 0.1+--+-- Default approximation ratio.+epsilon :: Double+epsilon = 0.1++-- | > (size * (1 - epsilon), size * (1 + epsilon))+tolerance :: Double -> Int -> (Int, Int)+tolerance epsilon size = (size - delta, size + delta)+  where+    delta = ceiling (fromIntegral size * epsilon)++-- * Auxiliary definitions++memo+  :: (t -> [t2] -> SG r)+  -> (SG r -> t1 -> Maybe Int -> a)+  -> t -> t1 -> Int -> a+memo make apply aliases k = generators+  where+    sg = make aliases []+    generators = sparseSized (apply sg k . Just) (99 <$ maxSizeM sg)++-- Oracles are computed only for sizes that are a power of two away from+-- the minimum size of the datatype @minSize + 2 ^ e@.+sparseSized :: (Int -> a) -> Maybe Int -> Int -> a+sparseSized f maxSizeM =+  maybe a0 snd . \size' -> find ((>= size') . fst) as+  where+    as = [ (s, f s) | s <- ss ]+    ss = 0 : maybe id (takeWhile . (>)) maxSizeM [ 2 ^ e | e <- [ 0 :: Int ..] ]+    a0 = f (fromJust maxSizeM)

@@ -0,0 +1,506 @@+-- | Internal module++{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE DeriveGeneric #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE ScopedTypeVariables #-}++module Boltzmann.Data.Oracle where++import Control.Applicative+import Control.Monad+import Control.Monad.Fix+import Control.Monad.Reader+import Control.Monad.State+import Data.Bifunctor+import Data.Data+import Data.Hashable ( Hashable )+import Data.HashMap.Lazy ( HashMap )+import qualified Data.HashMap.Lazy as HashMap+import Data.Maybe ( fromJust, isJust )+import Data.Monoid+import qualified Data.Vector as V+import GHC.Generics ( Generic )+import Numeric.AD+import Boltzmann.Data.Common+import Boltzmann.Data.Types+import Boltzmann.Solver++-- | We build a dictionary which reifies type information in order to+-- create a Boltzmann generator.+--+-- We denote by @n@ (or 'count') the number of types in the dictionary.+--+-- Every type has an index @0 <= i < n@; the variable @X i@ represents its+-- generating function @C_i(x)@, and @X (i + k*n)@ the GF of its @k@-th+-- "pointing" @C_i[k](x)@; we have+--+-- @+--   C_i[0](x) = C_i(x)+--   C_i[k+1](x) = x * C_i[k]'(x)+-- @+--+-- where @C_i[k]'@ is the derivative of @C_i[k]@. See also 'point'.+--+-- The /order/ (or /valuation/) of a power series is the index of the first+-- non-zero coefficient, called the /leading coefficient/.++data DataDef m = DataDef+  { count :: Int -- ^ Number of registered types+  , points :: Int -- ^ Number of iterations of the pointing operator+  , index :: HashMap TypeRep (Either Aliased Ix) -- ^ Map from types to indices+  , xedni :: HashMap Ix SomeData' -- ^ Inverse map from indices to types+  , xedni' :: HashMap Aliased (Ix, Alias m) -- ^ Inverse map to aliases+  , types :: HashMap C [(Integer, Constr, [C'])]+  -- ^ Structure of types and their pointings (up to 'points', initially 0)+  --+  -- Primitive types and empty types are mapped to an empty constructor list, and+  -- can be distinguished using 'Data.Data.dataTypeRep' on the 'SomeData'+  -- associated to it by 'xedni'.+  --+  -- The integer is a multiplicity which can be > 1 for pointings.+  , lTerm :: HashMap Ix (Nat, Integer)+  -- ^ Leading term @a * x ^ u@ of the generating functions @C_i[k](x)@ in the+  -- form (u, a).+  --+  -- [Order @u@] Smallest size of objects of a given type.+  -- [Leading coefficient @a@] number of objects of smallest size.+  , degree :: HashMap Ix Int+  -- ^ Degrees of the generating functions, when applicable: greatest size of+  -- objects of a given type.+  } deriving Show++-- | A pair @C i k@ represents the @k@-th "pointing" of the type at index @i@,+-- with generating function @C_i[k](x)@.+data C = C Ix Int+  deriving (Eq, Ord, Show, Generic)++instance Hashable C++data AC = AC Aliased Int+  deriving (Eq, Ord, Show, Generic)++instance Hashable AC++type C' = (Maybe Aliased, C)++newtype Aliased = Aliased Int+  deriving (Eq, Ord, Show, Generic)++instance Hashable Aliased++type Ix = Int++data Nat = Zero | Succ Nat+  deriving (Eq, Ord, Show)++instance Monoid Nat where+  mempty = Zero+  mappend (Succ n) = Succ . mappend n+  mappend Zero = id++natToInt :: Nat -> Int+natToInt Zero = 0+natToInt (Succ n) = 1 + natToInt n++infinity :: Nat+infinity = Succ infinity++dataDef :: [Alias m] -> DataDef m+dataDef as = DataDef+  { count = 0+  , points = 0+  , index = index+  , xedni = HashMap.empty+  , xedni' = xedni'+  , types = HashMap.empty+  , lTerm = HashMap.empty+  , degree = HashMap.empty+  } where+    xedni' = HashMap.fromList (fmap (\(i, a) -> (i, (-1, a))) as')+    index = HashMap.fromList (fmap (\(i, a) -> (ofType a, Left i)) as')+    as' = zip (fmap Aliased [0 ..]) as+    ofType (Alias f) = typeRep (f undefined)++-- | Find all types that may be types of subterms of a value of type @a@.+--+-- This will loop if there are infinitely many such types.+collectTypes :: Data a => [Alias m] -> proxy a -> DataDef m+collectTypes as a = collectTypesM a `execState` dataDef as++-- | Primitive datatypes have @C(x) = x@: they are considered as+-- having a single object (@lCoef@) of size 1 (@order@)).+primOrder :: Int+primOrder = 1++primOrder' :: Nat+primOrder' = Succ Zero++primlCoef :: Integer+primlCoef = 1++-- | The type of the first argument of 'Data.Data.gunfold'.+type GUnfold m = forall b r. Data b => m (b -> r) -> m r++-- | Type of 'xedni''.+type AMap m = HashMap Aliased (Ix, Alias m)++collectTypesM :: Data a => proxy a+  -> State (DataDef m) (Either Aliased Ix, ((Nat, Integer), Maybe Int))+collectTypesM a = chaseType a (const id)++chaseType :: Data a => proxy a+  -> ((Maybe (Alias m), Ix) -> AMap m -> AMap m)+  -> State (DataDef m) (Either Aliased Ix, ((Nat, Integer), Maybe Int))+chaseType a k = do+  let t = typeRep a+  dd@DataDef{..} <- get+  let+    lookup i r =+      let+        lTerm_i = lTerm #! i+        degree_i = HashMap.lookup i degree+      in return (r, (lTerm_i, degree_i))+  case HashMap.lookup t index of+    Nothing -> do+      let i = count+      put dd+        { count = i + 1+        , index = HashMap.insert t (Right i) index+        , xedni = HashMap.insert i (someData' a) xedni+        , xedni' = k (Nothing, i) xedni'+        }+      traverseType a i -- Updates lTerm and degree+    Just (Right i) -> do+      put dd { xedni' = k (Nothing, i) xedni' }+      lookup i (Right i)+    Just (Left j) ->+      case xedni' #! j of+        (-1, Alias f) -> do+          (_, ld) <- chaseType (ofType f) $ \(alias, i) ->+            let+              alias' = case alias of+                Nothing -> Alias f+                Just (Alias g) -> Alias (composeCastM f g)+            in+            k (Just alias', i) . HashMap.insert j (i, alias')+          return (Left j, ld)+        (i, _) -> lookup i (Left j)+  where+    ofType :: (m a -> m b) -> m a+    ofType _ = undefined++-- | Traversal of the definition of a datatype.+traverseType+  :: Data a => proxy a -> Ix+  -> State (DataDef m) (Either Aliased Ix, ((Nat, Integer), Maybe Int))+traverseType a i = do+  let d = withProxy dataTypeOf a+  mfix $ \ ~(_, (lTerm_i0, _)) -> do+    modify $ \dd@DataDef{..} -> dd+      { lTerm = HashMap.insert i lTerm_i0 lTerm+      }+    (types_i, ld@(_, degree_i)) <- traverseType' a d+    modify $ \dd@DataDef{..} -> dd+      { types = HashMap.insert (C i 0) types_i types+      , degree = maybe id (HashMap.insert i) degree_i degree+      }+    return (Right i, ld)++traverseType'+  :: Data a => proxy a -> DataType+  -> State (DataDef m)+      ([(Integer, Constr, [(Maybe Aliased, C)])], ((Nat, Integer), Maybe Int))+traverseType' a d | isAlgType d = do+  let+    constrs = dataTypeConstrs d+    collect+      :: GUnfold (StateT+        ([Either Aliased Ix], (Nat, Integer), Maybe Int)+        (State (DataDef m)))+    collect mkCon = do+      f <- mkCon+      let ofType :: (b -> a) -> Proxy b+          ofType _ = Proxy+          b = ofType f+      (j, (lTerm_, degree_)) <- lift (collectTypesM b)+      modify $ \(js, lTerm', degree') ->+        (j : js, lMul lTerm_ lTerm', liftA2 (+) degree_ degree')+      return (withProxy f b)+  tlds <- forM constrs $ \constr -> do+    (js, lTerm', degree') <-+      gunfold collect return constr `proxyType` a+        `execStateT` ([], (Zero, 1), Just 1)+    dd <- get+    let+      c (Left j) = (Just j, C (fst (xedni' dd #! j)) 0)+      c (Right i) = (Nothing, C i 0)+    return ((1, constr, [ c j | j <- js]), lTerm', degree')+  let+    (types_i, ls, ds) = unzip3 tlds+    lTerm_i = first Succ (lSum ls)+    degree_i = maxDegree ds+  return (types_i, (lTerm_i, degree_i))+traverseType' _ _ =+  return ([], ((primOrder', primlCoef), Just primOrder))++-- | If @(u, a)@ represents a power series of leading term @a * x ^ u@, and+-- similarly for @(u', a')@, this finds the leading term of their sum.+--+-- The comparison of 'Nat' is unrolled here for maximum laziness.+lPlus :: (Nat, Integer) -> (Nat, Integer) -> (Nat, Integer)+lPlus (Zero, lCoef) (Zero, lCoef') = (Zero, lCoef + lCoef')+lPlus (Zero, lCoef) _ = (Zero, lCoef)+lPlus _ (Zero, lCoef') = (Zero, lCoef')+lPlus (Succ order, lCoef) (Succ order', lCoef') =+  first Succ $ lPlus (order, lCoef) (order', lCoef')++-- | Sum of a list of series.+lSum :: [(Nat, Integer)] -> (Nat, Integer)+lSum [] = (infinity, 0)+lSum ls = foldl1 lPlus ls++-- | Leading term of a product of series.+lMul :: (Nat, Integer) -> (Nat, Integer) -> (Nat, Integer)+lMul (order, lCoef) (order', lCoef') = (order <> order', lCoef * lCoef')++lProd :: [(Nat, Integer)] -> (Nat, Integer)+lProd = foldl lMul (Zero, 1)++maxDegree :: [Maybe Int] -> Maybe Int+maxDegree = foldl (liftA2 max) (Just minBound)++-- | Pointing operator.+--+-- Populates a 'DataDef' with one more level of pointings.+-- ('collectTypes' produces a dictionary at level 0.)+--+-- The "pointing" of a type @t@ is a derived type whose values are essentially+-- values of type @t@, with one of their constructors being "pointed".+-- Alternatively, we may turn every constructor into variants that indicate+-- the position of points.+--+-- @+--   -- Original type+--   data Tree = Node Tree Tree | Leaf+--   -- Pointing of Tree+--   data Tree'+--     = Tree' Tree -- Point at the root+--     | Node'0 Tree' Tree -- Point to the left+--     | Node'1 Tree Tree' -- Point to the right+--   -- Pointing of the pointing+--   -- Notice that the "points" introduced by both applications of pointing+--   -- are considered different: exchanging their positions (when different)+--   -- produces a different tree.+--   data Tree''+--     = Tree'' Tree' -- Point 2 at the root, the inner Tree' places point 1+--     | Node'0' Tree' Tree -- Point 1 at the root, point 2 to the left+--     | Node'1' Tree Tree' -- Point 1 at the root, point 2 to the right+--     | Node'0'0 Tree'' Tree -- Points 1 and 2 to the left+--     | Node'0'1 Tree' Tree' -- Point 1 to the left, point 2 to the right+--     | Node'1'0 Tree' Tree' -- Point 1 to the right, point 2 to the left+--     | Node'0'1 Tree Tree'' -- Points 1 and 2 to the right+-- @+--+-- If we ignore points, some constructors are equivalent. Thus we may simply+-- calculate their multiplicity instead of duplicating them.+--+-- Given a constructor with @c@ arguments @C x_1 ... x_c@, and a sequence+-- @p_0 + p_1 + ... + p_c = k@ corresponding to a distribution of @k@ points+-- (@p_0@ are assigned to the constructor @C@ itself, and for @i > 0@, @p_i@+-- points are assigned within the @i@-th subterm), the multiplicity of the+-- constructor paired with that distribution is the multinomial coefficient+-- @multinomial k [p_1, ..., p_c]@.++point :: DataDef m -> DataDef m+point dd@DataDef{..} = dd+  { points = points'+  , types = foldl g types [0 .. count-1]+  } where+    points' = points + 1+    g types i = HashMap.insert (C i points') (types' i) types+    types' i = types #! C i 0 >>= h+    h (_, constr, js) = do+      ps <- partitions points' (length js)+      let+        mult = multinomial points' ps+        js' = zipWith (\(j', C i _) p -> (j', C i p)) js ps+      return (mult, constr, js')++-- | An oracle gives the values of the generating functions at some @x@.+type Oracle = HashMap C Double++-- | Find the value of @x@ such that the average size of the generator+-- for the @k-1@-th pointing is equal to @size@, and produce the associated+-- oracle. If the size is @Nothing@, find the radius of convergence.+--+-- The search evaluates the generating functions for some values of @x@ in+-- order to run a binary search. The evaluator is implemented using Newton's+-- method, the convergence of which has been shown for relevant systems in+-- /Boltzmann Oracle for Combinatorial Systems/,+-- C. Pivoteau, B. Salvy, M. Soria.+makeOracle :: DataDef m -> TypeRep -> Maybe Double -> Oracle+makeOracle dd0 t size' =+  seq v+  HashMap.fromList (zip cs (V.toList v))+  where+    -- We need the next pointing to capture the average size in an equation.+    dd@DataDef{..} = if isJust size' then point dd0 else dd0+    cs = flip C <$> [0 .. points] <*> [0 .. count - 1]+    m = count * (points + 1)+    k = points - 1+    i = case index #! t of+      Left j -> fst (xedni' #! j)+      Right i -> i+    checkSize _ (Just ys) | V.any (< 0) ys = False+    -- There may be solutions outside of the radius+    -- of convergence, but with negative components.+    checkSize (Just size) (Just ys) =+      size >= size_+      where+        size_ = ys V.! j' / ys V.! j+        j = dd ? C i k+        j' = dd ? C i (k + 1)+    checkSize Nothing (Just _) = True+    checkSize _ Nothing = False+    -- Equations defining C_i(x) for all types with indices i+    phis :: Num a => V.Vector (a -> V.Vector a -> a)+    phis = V.fromList [ phi dd c (types #! c) | c <- listCs dd ]+    eval' :: Double -> Maybe (V.Vector Double)+    eval' x = fixedPoint defSolveArgs phi' (V.replicate m 0)+      where+        phi' :: (Mode a, Scalar a ~ Double) => V.Vector a -> V.Vector a+        phi' y = fmap (\f -> f (auto x) y) phis+    v = (fromJust . snd) (search eval' (checkSize size'))++-- | Generating function definition. This defines a @Phi_i[k]@ function+-- associated with the @k@-th pointing of the type at index @i@, such that:+--+-- > C_i[k](x)+-- >   = Phi_i[k](x, C_0[0](x), ..., C_(n-1)[0](x),+-- >              ..., C_0[k](x), ..., C_(n-1)[k](x))+--+-- Primitive datatypes have @C(x) = x@: they are considered as+-- having a single object ('lCoef') of size 1 ('order')).+phi :: Num a => DataDef m -> C -> [(Integer, constr, [C'])]+  -> a -> V.Vector a -> a+phi DataDef{..} (C i _) [] =+  case xedni #! i of+    SomeData a ->+      case (dataTypeRep . withProxy dataTypeOf) a of+        AlgRep _ -> \_ _ -> 0+        _ -> \x _ -> fromInteger primlCoef * x ^ primOrder+phi dd@DataDef{..} _ tyInfo = f+  where+    f x y = x * (sum . fmap (toProd y)) tyInfo+    toProd y (w, _, js) =+      fromInteger w * product [ y V.! (dd ? j) | (_, j) <- js ]++-- | Maps a key representing a type @a@ (or one of its pointings) to a+-- generator @m a@.+type Generators m = (HashMap AC (SomeData m), HashMap C (SomeData m))++-- | Build all involved generators at once.+makeGenerators+  :: forall m. MonadRandomLike m+  => DataDef m -> Oracle -> Generators m+makeGenerators DataDef{..} oracle =+  seq oracle+  (generatorsL, generatorsR)+  where+    f (C i _) tyInfo = case xedni #! i of+      SomeData a -> SomeData $ incr >>+        case tyInfo of+          [] -> defGen+          _ -> frequencyWith doubleR (fmap g tyInfo) `proxyType` a+    g :: Data a => (Integer, Constr, [C']) -> (Double, m a)+    g (v, constr, js) =+      ( fromInteger v * w+      , gunfold generate return constr `runReaderT` gs)+      where+        gs = fmap (\(j', i) -> m j' i) js+        m = maybe (generatorsR #!) m'+        m' j (C _ k) = (generatorsL #! AC j k)+        w = product $ fmap ((oracle #!) . snd) js+    h (j, (i, Alias f)) k =+      (AC j k, applyCast f (generatorsR #! C i k))+    generatorsL = HashMap.fromList (liftA2 h (HashMap.toList xedni') [0 .. points])+    generatorsR = HashMap.mapWithKey f types++type SmallGenerators m =+  (HashMap Aliased (SomeData m), HashMap Ix (SomeData m))++-- | Generators of values of minimal sizes.+smallGenerators+  :: forall m. MonadRandomLike m => DataDef m -> SmallGenerators m+smallGenerators DataDef{..} = (generatorsL, generatorsR)+  where+    f i (SomeData a) = SomeData $ incr >>+      case types #! C i 0 of+        [] -> defGen+        tyInfo ->+          let gs = (tyInfo >>= g (fst (lTerm #! i))) in+          frequencyWith integerR gs `proxyType` a+    g :: Data a => Nat -> (Integer, Constr, [C']) -> [(Integer, m a)]+    g minSize (_, constr, js) =+      guard (minSize == Succ size) *>+      [(weight, gunfold generate return constr `runReaderT` gs)]+      where+        (size, weight) = lProd [ lTerm #! i | (_, C i _) <- js ]+        gs = fmap lookup js+        lookup (j', C i _) = maybe (generatorsR #! i) (generatorsL #!) j'+    h (j, (i, Alias f)) = (j, applyCast f (generatorsR #! i))+    generatorsL = (HashMap.fromList . fmap h . HashMap.toList) xedni'+    generatorsR = HashMap.mapWithKey f xedni++generate :: Applicative m => GUnfold (ReaderT [SomeData m] m)+generate rest = ReaderT $ \(g : gs) ->+  rest `runReaderT` gs <*> unSomeData g++defGen :: (Data a, MonadRandomLike m) => m a+defGen = gen+  where+    gen =+      let dt = withProxy dataTypeOf gen in+      case dataTypeRep dt of+        IntRep -> fromConstr . mkIntegralConstr dt <$> int+        FloatRep -> fromConstr . mkRealConstr dt <$> double+        CharRep -> fromConstr . mkCharConstr dt <$> char+        AlgRep _ -> error "Cannot generate for empty type."+        NoRep -> error "No representation."++-- * Short operators++(?) :: DataDef m -> C -> Int+dd ? C i k = i + k * count dd++-- | > dd ? (listCs dd !! i) = i+listCs :: DataDef m -> [C]+listCs dd = liftA2 (flip C) [0 .. points dd] [0 .. count dd - 1]++ix :: C -> Int+ix (C i _) = i++-- | > dd ? (dd ?! i) = i+(?!) :: DataDef m -> Int -> C+dd ?! j = C i k+  where (k, i) = j `divMod` count dd++getGenerator :: Data a => DataDef m -> Generators m -> proxy a -> Int -> m a+getGenerator dd (l, r) a k = unSomeData $+  case index dd #! typeRep a of+    Right i -> (r #! C i k)+    Left j -> (l #! AC j k)++getSmallGenerator :: Data a => DataDef m -> SmallGenerators m -> proxy a -> m a+getSmallGenerator dd (l, r) a = unSomeData $+  case index dd #! typeRep a of+    Right i -> (r #! i)+    Left j -> (l #! j)++(#!) :: (Eq k, Hashable k)+  => HashMap k v -> k -> v+(#!) = (HashMap.!)

@@ -0,0 +1,197 @@+-- | Internal module++{-# LANGUAGE GADTs #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeOperators #-}++module Boltzmann.Data.Types where++import Control.Monad.Random+import Control.Monad.Trans+import Data.Coerce+import Data.Data+import Data.Function+import Test.QuickCheck++data SomeData m where+  SomeData :: Data a => m a -> SomeData m++type SomeData' = SomeData Proxy++-- | Dummy instance for debugging.+instance Show (SomeData m) where+  show _ = "SomeData"++data Alias m where+  Alias :: (Data a, Data b) => !(m a -> m b) -> Alias m++type AliasR m = Alias (RejectT m)++-- | Dummy instance for debugging.+instance Show (Alias m) where+  show _ = "Alias"++-- | Main constructor for 'Alias'.+alias :: (Monad m, Data a, Data b) => (a -> m b) -> Alias m+alias = Alias . (=<<)++-- | Main constructor for 'AliasR'.+aliasR :: (Monad m, Data a, Data b) => (a -> m b) -> AliasR m+aliasR = Alias . (=<<) . fmap lift++-- | > coerceAlias :: Alias m -> Alias (AMonadRandom m)+coerceAlias :: Coercible m n => Alias m -> Alias n+coerceAlias = coerce++-- | > coerceAliases :: [Alias m] -> [Alias (AMonadRandom m)]+coerceAliases :: Coercible m n => [Alias m] -> [Alias n]+coerceAliases = coerce++-- | > composeCast f g = f . g+composeCastM :: forall a b c d m+  . (Typeable b, Typeable c)+  => (m c -> d) -> (a -> m b) -> (a -> d)+composeCastM f g | Just Refl <- eqT :: Maybe (b :~: c) = f . g+composeCastM _ _ = castError ([] :: [b]) ([] :: [c])++castM :: forall a b m+  . (Typeable a, Typeable b)+  => m a -> m b+castM a | Just Refl <- eqT :: Maybe (a :~: b) = a+castM a = let x = castError a x in x++unSomeData :: Typeable a => SomeData m -> m a+unSomeData (SomeData a) = castM a++applyCast :: (Typeable a, Data b) => (m a -> m b) -> SomeData m -> SomeData m+applyCast f = SomeData . f . unSomeData++castError :: (Typeable a, Typeable b)+  => proxy a -> proxy' b -> c+castError a b = error $ unlines+  [ "Error trying to cast"+  , "  " ++ show (typeRep a)+  , "to"+  , "  " ++ show (typeRep b)+  ]++withProxy :: (a -> b) -> proxy a -> b+withProxy f _ =+  f (error "This should not be evaluated\n")++reproxy :: proxy a -> Proxy a+reproxy _ = Proxy++proxyType :: m a -> proxy a -> m a+proxyType = const++someData' :: Data a => proxy a -> SomeData'+someData' = SomeData . reproxy++-- | Size as the number of constructors.+type Size = Int++-- | Internal transformer for rejection sampling.+--+-- > ReaderT Size (StateT Size (MaybeT m)) a+newtype RejectT m a = RejectT+  { unRejectT :: forall r. Size -> Size -> m r -> (Size -> a -> m r) -> m r+  }++instance Functor (RejectT m) where+  fmap f (RejectT go) = RejectT $ \maxSize size retry cont ->+    go maxSize size retry $ \size a -> cont size (f a)++instance Applicative (RejectT m) where+  pure a = RejectT $ \_maxSize size _retry cont ->+    cont size a+  RejectT f <*> RejectT x = RejectT $ \maxSize size retry cont ->+    f maxSize size retry $ \size f_ ->+      x maxSize size retry $ \size x_ ->+        cont size (f_ x_)++instance Monad (RejectT m) where+  RejectT x >>= f = RejectT $ \maxSize size retry cont ->+    x maxSize size retry $ \size x_ ->+      unRejectT (f x_) maxSize size retry cont++instance MonadTrans RejectT where+  lift m = RejectT $ \_maxSize size _retry cont ->+    m >>= cont size++-- | Set lower bound+runRejectT :: Monad m => (Size, Size) -> RejectT m a -> m a+runRejectT (minSize, maxSize) (RejectT m) = fix $ \go ->+  m maxSize 0 go $ \size a ->+    if size < minSize then+      go+    else+      return a+--runRejectT (minSize, maxSize) (RejectT m) = fix $ \go -> do+--  x' <- runMaybeT (m `runReaderT` maxSize `runStateT` 0)+--  case x' of+--    Just (x, size) | size >= minSize -> return x+--    _ -> go++newtype AMonadRandom m a = AMonadRandom+  { asMonadRandom :: m a+  } deriving (Functor, Applicative, Monad)++instance MonadTrans AMonadRandom where+  lift = AMonadRandom++-- ** Dictionaries++-- | @'MonadRandomLike' m@ defines basic components to build generators,+-- allowing the implementation to remain abstract over both the+-- 'Test.QuickCheck.Gen' type and 'MonadRandom' instances.+--+-- For the latter, the wrapper 'AMonadRandom' is provided to avoid+-- overlapping instances.+class Monad m => MonadRandomLike m where+  -- | Called for every constructor. Counter for ceiled rejection sampling.+  incr :: m ()+  incr = return ()++  -- | @doubleR upperBound@: generates values in @[0, upperBound]@.+  doubleR :: Double -> m Double++  -- | @integerR upperBound@: generates values in @[0, upperBound-1]@.+  integerR :: Integer -> m Integer++  -- | Default @Int@ generator.+  int :: m Int++  -- | Default @Double@ generator.+  double :: m Double++  -- | Default @Char@ generator.+  char :: m Char++instance MonadRandomLike Gen where+  doubleR x = choose (0, x)+  integerR x = choose (0, x-1)+  int = arbitrary+  double = arbitrary+  char = arbitrary++instance MonadRandomLike m => MonadRandomLike (RejectT m) where+  incr = RejectT $ \maxSize size retry cont ->+    if size >= maxSize then+      retry+    else+      cont (size + 1) ()+  doubleR = lift . doubleR+  integerR = lift . integerR+  int = lift int+  double = lift double+  char = lift char++instance MonadRandom m => MonadRandomLike (AMonadRandom m) where+  doubleR x = lift $ getRandomR (0, x)+  integerR x = lift $ getRandomR (0, x-1)+  int = lift getRandom+  double = lift getRandom+  char = lift getRandom

@@ -0,0 +1,69 @@+-- | Solve systems of equations++{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE RankNTypes #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE TypeFamilies #-}++module Boltzmann.Solver where++import Control.Applicative+import Data.AEq ( (~==) )+import Numeric.AD.Mode+import Numeric.AD.Mode.Forward+import Numeric.LinearAlgebra+import qualified Data.Vector as V+import qualified Data.Vector.Storable as S++data SolveArgs = SolveArgs+  { accuracy :: Double+  , numIterations :: Int+  } deriving (Eq, Ord, Show)++defSolveArgs :: SolveArgs+defSolveArgs = SolveArgs 1e-8 20++findZero+  :: SolveArgs+  -> (forall s. V.Vector (AD s (Forward R)) -> V.Vector (AD s (Forward R)))+  -> Vector R+  -> Maybe (Vector R)+findZero SolveArgs{..} f = newton numIterations+  where+    newton 0 _ = Nothing+    newton n x+      | norm_y == 1/0 = Nothing+      | norm_y > accuracy = newton (n - 1) (x - jacobian <\> y)+      | otherwise = Just x+      where+        norm_y = norm_Inf y+        jacobian = (fromRows . V.toList . fmap (V.convert . snd)) yj+        y = (V.convert . fmap fst) yj+        yj = jacobian' f (S.convert x)++fixedPoint+  :: SolveArgs+  -> (forall a. (Mode a, Scalar a ~ R) => V.Vector a -> V.Vector a)+  -> V.Vector R+  -> Maybe (V.Vector R)+fixedPoint args f =+  fmap S.convert . findZero args (liftA2 (V.zipWith (-)) f id) . S.convert++-- | Assuming @p . f@ is satisfied only for positive values in some interval+-- @(0, r]@, find @f r@.+search :: (Double -> a) -> (a -> Bool) -> (Double, a)+search f p = search' e0 (0 : [2 ^ n | n <- [0 .. 100 :: Int]])+  where+    search' y (x : xs@(x' : _))+      | p y' = search' y' xs+      | otherwise = search'' y x x'+      where y' = f x'+    search' _ _ = error "Solution not found. Uncontradictable predicate?"+    search'' y x x'+      | x ~== x' = (x, y)+      | p y_ = search'' y_ x_ x'+      | otherwise = search'' y x x_+      where+        x_ = (x + x') / 2+        y_ = f x_+    e0 = error "Solution not found. Unsatisfiable predicate?"

@@ -0,0 +1,220 @@+-- | Applicative interface to define recursive structures and derive Boltzmann+-- samplers.+--+-- Given the recursive structure of the types, and how to combine generators,+-- the library takes care of computing the oracles and setting the right+-- distributions.++{-# LANGUAGE FlexibleContexts, FlexibleInstances, GADTs, RankNTypes #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE DeriveFunctor, DeriveGeneric, ImplicitParams #-}+{-# LANGUAGE RecordWildCards, DeriveDataTypeable #-}+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses #-}++module Boltzmann.Species where++import Control.Applicative+import Control.Monad+import Data.Bifunctor+import Data.Coerce+import Data.Function+import Data.Foldable+import Data.List+import Data.Maybe+import Data.Vector ( Vector )+import qualified Data.Vector as V+import qualified Numeric.AD as AD++import Boltzmann.Data.Common+import Boltzmann.Data.Types+import Boltzmann.Solver++class Embed f m where+  emap :: (m a -> m b) -> f a -> f b+  -- | A natural transformation between @f@ and @m@?+  embed :: m a -> f a++-- | 'Applicative' defines a product, 'Alternative' defines an addition,+-- with scalar multiplication we get a module.+--+-- This typeclass allows to directly tweak weights in the oracle by+-- chosen factors.+class (Alternative f, Num (Scalar f)) => Module f where+  type Scalar f :: *++  -- | Scalar embedding.+  scalar :: Scalar f -> f ()+  scalar x = x <.> pure ()++  -- | Scalar multiplication.+  (<.>) :: Scalar f -> f a -> f a+  x <.> f = scalar x *> f++infixr 3 <.>++type Endo a = a -> a++data System f a c = System+  { dim :: Int+  , sys' :: f () -> Vector (f a) -> (Vector (f a), c)+  } deriving (Functor)++sys :: System f a c -> f () -> Vector (f a) -> Vector (f a)+sys = (fmap . fmap . fmap) fst sys'++newtype ConstModule r a = ConstModule { unConstModule :: r }++instance Functor (ConstModule r) where+  fmap _ (ConstModule r) = ConstModule r++instance Num r => Embed (ConstModule r) m where+  emap _ (ConstModule r) = ConstModule r+  embed _ = ConstModule 1++instance Num r => Applicative (ConstModule r) where+  pure _ = ConstModule 1+  ConstModule x <*> ConstModule y = ConstModule (x * y)++instance Num r => Alternative (ConstModule r) where+  empty = ConstModule 0+  ConstModule x <|> ConstModule y = ConstModule (x + y)++instance Num r => Module (ConstModule r) where+  type Scalar (ConstModule r) = r+  scalar = ConstModule+  x <.> ConstModule r = ConstModule (x * r)++solve+  :: forall b c+  . (forall a. Num a => System (ConstModule a) b c)+  -> Double -> Maybe (Vector Double)+solve s x = fixedPoint defSolveArgs phi' (V.replicate (dim s') 0)+  where+    phi' :: forall a. (AD.Mode a, AD.Scalar a ~ Double) => Endo (Vector a)+    phi' = coerce (sys s (scalar (AD.auto x)) :: Endo (Vector (ConstModule a b)))+    -- Arbitrary instantiation to get its dimension.+    s' :: System (ConstModule Int) b c+    s' = s++sizedGenerator+  :: forall b c m+  . MonadRandomLike m+  => (forall f. (Module f, Embed f m) => System (Pointiful f) b c)+  -> Int  -- ^ Index of type+  -> Int  -- ^ Points+  -> Maybe Double  -- ^ Expected size (or singular sampler)+  -> m b+sizedGenerator s i k size' = fst (sfix s' x oracle) V.! j+  where+    (x, oracle) = solveSized s i k size'+    s' = point (k + 1) s+    j = i * (k + 2) + k++solveSized+  :: forall b c+  . (forall a. Num a => System (Pointiful (ConstModule a)) b c)+  -> Int  -- ^ Index of type+  -> Int  -- ^ Points+  -> Maybe Double  -- ^ Expected size (or singular sampler)+  -> (Double, Vector Double)+solveSized s i k size' =+  fmap fromJust (search (solve s') (checkSize size'))+  where+    s' :: forall a. Num a => System (ConstModule a) b c+    s' = point (k + 1) s+    j = i * (k + 2) + k+    j' = i * (k + 2) + k + 1+    checkSize _ (Just ys) | V.any (< 0) ys = False+    checkSize (Just size) (Just ys) = size >= ys V.! j' / ys V.! j+    checkSize Nothing (Just _) = True+    checkSize _ Nothing = False++newtype Weighted m a = Weighted [(Double, m a)]++weighted :: Double -> m a -> Weighted m a+weighted x a = Weighted [(x, a)]++runWeighted :: MonadRandomLike m => Weighted m a -> (Double, m a)+runWeighted (Weighted [a]) = a+runWeighted (Weighted as) = (sum (fmap fst as), frequencyWith doubleR as)++instance Functor m => Functor (Weighted m) where+  fmap f (Weighted as) = Weighted ((fmap . fmap . fmap) f as)++instance MonadRandomLike m => Embed (Weighted m) m where+  emap f = Weighted . (: []) . fmap f . runWeighted+  embed m = Weighted [(1, m)]++instance MonadRandomLike m => Applicative (Weighted m) where+  pure a = Weighted [(1, pure a)]+  f' <*> a' = Weighted [(u * v, f <*> a)]+    where+      (u, f) = runWeighted f'+      (v, a) = runWeighted a'++instance MonadRandomLike m => Alternative (Weighted m) where+  empty = Weighted []+  Weighted as <|> Weighted bs = Weighted (as ++ bs)++instance MonadRandomLike m => Module (Weighted m) where+  type Scalar (Weighted m) = Double+  scalar x = Weighted [(x, pure ())]+  x <.> Weighted as = Weighted (fmap (first (x *)) as)++sfix+  :: MonadRandomLike m+  => System (Weighted m) b c -> Double -> Vector Double -> (Vector (m b), c)+sfix s x oracle =+  fix $+    (first . fmap) (snd . runWeighted) .+    sys' s (scalar x) .+    V.zipWith weighted oracle .+    fst++data Pointiful f a = Pointiful [f a] | Zero (f a)++instance Functor f => Functor (Pointiful f) where+  fmap f (Pointiful v) = Pointiful ((fmap . fmap) f v)+  fmap f (Zero x) = Zero (fmap f x)++instance Embed f m => Embed (Pointiful f) m where+  emap f (Pointiful v) = Pointiful ((fmap . emap) f v)+  emap f (Zero x) = Zero (emap f x)+  embed = Zero . embed++instance Module f => Applicative (Pointiful f) where+  pure a = Zero (pure a)+  Zero f <*> Zero x = Zero (f <*> x)+  Zero f <*> Pointiful xs = Pointiful (fmap (f <*>) xs)+  Pointiful fs <*> Zero x = Pointiful (fmap (<*> x) fs)+  Pointiful fs <*> Pointiful xs = Pointiful (convolute fs xs)+    where+      convolute fs xs = zipWith3 sumOfProducts [0 ..] (inits' fs) (inits' xs)+      inits' = tail . inits+      sumOfProducts k f x = asum (zipWith3 (times k) [0 ..] f (reverse x))+      times k k1 f x = fromInteger (binomial k k1) <.> f <*> x++instance Module f => Alternative (Pointiful f) where+  empty = Zero empty+  Pointiful xs <|> Pointiful ys = Pointiful (zipWith (<|>) xs ys)+  Pointiful (x : xs) <|> Zero y = Pointiful ((x <|> y) : xs)+  Zero x <|> Pointiful (y : ys) = Pointiful ((x <|> y) : ys)+  Zero x <|> Zero y = Zero (x <|> y)+  Pointiful [] <|> m = m+  m <|> Pointiful [] = m++instance Module f => Module (Pointiful f) where+  type Scalar (Pointiful f) = Scalar f+  scalar = Zero . scalar++unPointiful :: Alternative f => Pointiful f a -> [f a]+unPointiful (Pointiful as) = as+unPointiful (Zero a) = a : repeat empty++point :: Module f => Int -> System (Pointiful f) b c -> System f b c+point k s = System ((k + 1) * dim s) $ \x ->+  first flatten . sys' s (Pointiful (repeat x)) . resize+  where+    flatten = join . fmap (V.fromList . take (k + 1) . unPointiful)+    resize v = V.generate (dim s) $ \i ->+      Pointiful [v V.! j | j <- [i * (k + 1) .. i * (k + 1) + k]]

@@ -0,0 +1,76 @@+module Test.Stats where++import Control.Monad++import Data.List+import Data.Maybe++mean :: Foldable v => v Int -> Double+mean xs = fromIntegral (sum xs) / fromIntegral (length xs)++-- | Number of samples to estimate a probability distribution on a finite set+-- of size @n@ to precision @epsilon@ (infinity-norm between distributions)+-- with probability at least @(1 - delta)@.+sampleSize+  :: Int  -- ^ Domain size+  -> Double  -- ^ Target distance (infinity-norm)+  -> Double  -- ^ Target error probability+  -> Int+sampleSize n epsilon delta =+  ceiling (log (2 * fromIntegral n / delta) / (2 * epsilon ^ 2))++-- | Number of trees with @n@ internal nodes.+catalan :: [Integer]+catalan = fmap catalan' [0 ..]+  where+    catalan' 0 = 1+    catalan' i =+      let prefix = take i catalan+      in sum $ zipWith (*) prefix (reverse prefix)++-- | Average size of a binary tree given the probability (@> 1/2@) of choosing+-- a leaf.+avgSize :: Fractional a => a -> a+avgSize p = 1 / (2 * p - 1)++-- | Inverse of 'avgSize'.+invAvgSize :: Fractional a => a -> a+invAvgSize s = (1 / s + 1) / 2++-- | Distribution of sizes (actually, @(size - 1) / 2@), given the probability+-- of choosing a leaf.+distribution :: Fractional a => a -> [a]+distribution p = zipWith f [0 ..] catalan+  where+    f i c = fromInteger c * p * (p * (1 - p)) ^ i++expected :: Fractional a => Maybe a -> (Int, Int) -> Double -> Double -> (Int, [(Int, a)])+expected avgSize' (minSize_, maxSize_) epsilon delta = (k, d)+  where+    p = maybe (1/2) invAvgSize avgSize'+    minSize = (minSize_ + 1) `div` 2+    maxSize = maxSize_ `div` 2+    n = maxSize - minSize + 1+    k = sampleSize n epsilon delta+    d_ = (take n . drop minSize . distribution) p+    d = zip [minSize ..] (fmap (/ sum d_) d_)++runExperiment+  :: (Fractional a, Ord a, Monad m)+  => (Int, [(Int, a)]) -> m Int -> m ([(Int, a)], [(Int, a)], a)+runExperiment (k, d) gen = cmp' . collect <$> replicateM k gen+  where+    collect :: Fractional a => [Int] -> [(Int, a)]+    collect = fmap c . group . sort+    c xs@(x : _) = (x, fromIntegral (length xs) / fromIntegral k)+    c _ = undefined+    cmp' z = (d, z, cmp d z)+    cmp :: (Ord a, Num a) => [(Int, a)] -> [(Int, a)] -> a+    cmp xs ys = maximum (zipWith_ (\x y -> abs (x - y)) xs ys)+    zipWith_ :: (a -> a -> a) -> [(Int, a)] -> [(Int, a)] -> [a]+    zipWith_ f xxs@((x, m) : xs) yys@((y, n) : ys)+      | x == y = f m n : zipWith_ f xs ys+      | x < y = m : zipWith_ f xs yys+      | otherwise = n : zipWith_ f xxs ys+    zipWith_ f [] ys = fmap snd ys+    zipWith_ f xs [] = fmap snd xs

@@ -0,0 +1,85 @@+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE DeriveDataTypeable #-}+{-# LANGUAGE OverloadedStrings #-}+{-# LANGUAGE TypeOperators #-}++import Control.Monad+import Data.Data+import Data.Foldable+import Data.IORef+import Data.List+import System.Exit+import System.IO++import Options.Generic++import Boltzmann.Data+import Boltzmann.Data.Data++import Test.Stats++data T = L | N T T+  deriving (Eq, Ord, Show, Data)++size :: T -> Int+size (N l r) = 1 + size l + size r+size L = 0++eps, del :: Double+eps = 0.01+del = 0.001++-- | Periodically print stuff so that Travis does not think we're stuck.+counting x gen = do+  modifyIORef x (+ 1)+  readIORef x >>= \x ->+    when (x `mod` 1000 == 0) $ putStr "." >> hFlush stdout+  gen++-- | Invocation: stack test [--test-arguments TEST_SIZE]+type Input = Maybe (Int <?> "Test size")++main = do+  n_ <- getRecord "Test program" :: IO Input+  success <- newIORef True++  let n = maybe 10 unHelpful n_+      range = tolerance epsilon n++  for_+    [ ( "reject "+      , generatorSR+      , expected Nothing range eps del+      )+    , ( "rejectSimple "+      , generatorR'+      , expected (Just (fromIntegral n)) range eps del+      )+    ] $ \(name, g, kdist) -> do+    putStrLn $ name ++ show n+    let gen = (fmap size . asMonadRandom . g) n+    x <- newIORef 0+    (expectedDist, estimatedDist, diff) <- runExperiment kdist (counting x gen)+    putStrLn ""+    when (diff > eps) $ do+      writeIORef success False+      putStrLn $ "FAIL > " ++ show diff+      print expectedDist+      print estimatedDist++{-+  let k = 80000+      eps = 0.1+      gen = (fmap size . asMonadRandom . generatorP') n+  putStrLn $ "pointed " ++ show n+  x <- newIORef 0+  sizes <- replicateM k (counting x gen)+  putStrLn ""+  let diff = abs (mean sizes - fromIntegral (n `div` 2))+  when (diff > eps) $ do+    writeIORef success False+    putStrLn $ "FAIL > " ++ show diff+-}++  success <- readIORef success+  unless success exitFailure