diff --git a/README.md b/README.md --- a/README.md +++ b/README.md @@ -1,13 +1,17 @@ Generic random generators [![Hackage](https://img.shields.io/hackage/v/generic-random.svg)](https://hackage.haskell.org/package/generic-random) [![Build Status](https://travis-ci.org/Lysxia/generic-random.svg)](https://travis-ci.org/Lysxia/generic-random) ========================= +`Generic.Random.Data` +--------------------- + Define sized random generators for almost any type. ```haskell {-# LANGUAGE DeriveDataTypeable #-} + import Data.Data import Test.QuickCheck - import Data.Random.Generics + import Generic.Random.Data data Term = Lambda Int Term | App Term Term | Var Int deriving (Show, Data) @@ -25,11 +29,67 @@ - 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. -- Can be extended or modified with user defined generators. +- Uses `Data.Data` generics. +- Works with QuickCheck and MonadRandom, but also similar user-defined monads + for randomness (just implement `MonadRandomLike`). +- Can be tweaked somewhat with user defined generators. +`Generic.Random.Generic` +------------------------ + +Say goodbye to `Constructor <$> arbitrary <*> arbitrary <*> arbitrary`-boilerplate. + +```haskell + {-# LANGUAGE DataKinds #-} + {-# LANGUAGE DeriveGeneric #-} + {-# LANGUAGE TypeApplications #-} + + import GHC.Generics ( Generic ) + import Test.QuickCheck + import Generic.Random.Generic + + data Tree a = Leaf | Node (Tree a) a (Tree a) + deriving (Show, Generic) + + instance Arbitrary a => Arbitrary (Tree a) where + arbitrary = genericArbitrary' @'Z + + -- Equivalent to + -- > arbitrary = + -- > sized $ \n -> + -- > if n == 0 then + -- > return Leaf + -- > else + -- > oneof + -- > [ return Leaf + -- > , Node <$> arbitrary <*> arbitrary <*> arbitrary + -- > ] + + main = sample (arbitrary :: Gen (Tree ())) +``` + +- User-specified distribution of constructors. +- A simple (optional) strategy to ensure termination: `Test.QuickCheck.Gen`'s + size parameter decreases at every recursive `genericArbitrary'` call; when it + reaches zero, sample directly from a finite set of finite values. +- Uses `GHC.Generics` generics. +- Just for QuickCheck's `arbitrary`. +- More flexible than `Generic.Random.Data`'s Boltzmann samplers, which compute + fixed weights for a given target size and concrete type, but with a less + regular distribution. + +`Generic.Random.Boltzmann` +-------------------------- + +An experimental interface to obtain Boltzmann samplers from an applicative +specification of a combinatorial system. + +No documentation (yet). + References ---------- + +Papers about Boltzmann samplers, used in `Generic.Random.Data`: - 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), diff --git a/bench/binaryTree.hs b/bench/binaryTree.hs --- a/bench/binaryTree.hs +++ b/bench/binaryTree.hs @@ -14,9 +14,9 @@ import Test.QuickCheck.Gen import Test.QuickCheck.Random import Control.Exception ( evaluate ) -import Data.Random.Generics -import Data.Random.Generics.Internal -import Data.Random.Generics.Internal.Types +import Generic.Random.Data +import Generic.Random.Internal.Data +import Generic.Random.Internal.Types data T = N T T | L deriving (Eq, Ord, Show, Data, Generic) diff --git a/generic-random.cabal b/generic-random.cabal --- a/generic-random.cabal +++ b/generic-random.cabal @@ -1,7 +1,7 @@ name: generic-random -version: 0.1.1.0 +version: 0.2.0.0 synopsis: Generic random generators -description: Please see the README below. +description: Please see the README. homepage: http://github.com/lysxia/generic-random license: MIT license-file: LICENSE @@ -12,16 +12,20 @@ build-type: Simple extra-source-files: README.md cabal-version: >=1.10 -tested-with: GHC == 7.10.3 +tested-with: GHC == 8.0.1 library hs-source-dirs: src exposed-modules: - Data.Random.Generics - Data.Random.Generics.Internal - Data.Random.Generics.Internal.Oracle - Data.Random.Generics.Internal.Solver - Data.Random.Generics.Internal.Types + Generic.Random.Boltzmann + Generic.Random.Data + Generic.Random.Generic + Generic.Random.Internal.Common + Generic.Random.Internal.Data + Generic.Random.Internal.Generic + Generic.Random.Internal.Oracle + Generic.Random.Internal.Solver + Generic.Random.Internal.Types build-depends: base >= 4.8 && < 5, containers, diff --git a/src/Data/Random/Generics.hs b/src/Data/Random/Generics.hs deleted file mode 100644 --- a/src/Data/Random/Generics.hs +++ /dev/null @@ -1,302 +0,0 @@ --- | 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 Data.Random.Generics ( - 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 Data.Random.Generics.Internal -import Data.Random.Generics.Internal.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 @generic-random@ 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 get positive @Int@s: --- --- @ --- let --- as = ['alias' $ \\() -> 'choose' (0, 100) :: 'Gen' Int)] --- in --- 'generatorPWith' as 'asGen' :: 'Size' -> 'Gen' [Int] --- @ - -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 []) diff --git a/src/Data/Random/Generics/Internal.hs b/src/Data/Random/Generics/Internal.hs deleted file mode 100644 --- a/src/Data/Random/Generics/Internal.hs +++ /dev/null @@ -1,146 +0,0 @@ -{-# LANGUAGE RecordWildCards, DeriveFunctor #-} -module Data.Random.Generics.Internal where - -import Control.Arrow ( (&&&) ) -import Control.Applicative -import Data.Data -import Data.Foldable -import Data.Maybe -import qualified Data.HashMap.Lazy as HashMap -import Data.Random.Generics.Internal.Oracle -import Data.Random.Generics.Internal.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) diff --git a/src/Data/Random/Generics/Internal/Oracle.hs b/src/Data/Random/Generics/Internal/Oracle.hs deleted file mode 100644 --- a/src/Data/Random/Generics/Internal/Oracle.hs +++ /dev/null @@ -1,539 +0,0 @@ -{-# LANGUAGE FlexibleContexts, GADTs, RankNTypes, ScopedTypeVariables #-} -{-# LANGUAGE DeriveGeneric, ImplicitParams #-} -{-# LANGUAGE RecordWildCards, DeriveDataTypeable #-} -module Data.Random.Generics.Internal.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 qualified Data.Vector.Storable as S -import GHC.Generics ( Generic ) -import Numeric.AD -import Data.Random.Generics.Internal.Types -import Data.Random.Generics.Internal.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 (S.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) | S.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 S.! j' / ys S.! 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 (S.Vector Double) - eval' x = fixedPoint defSolveArgs phi' (S.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 (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 :: (Functor m, 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 :: (Functor m, 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) - --- * General helper functions - -frequencyWith - :: (Show r, 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. - -(#!) :: (Eq k, Hashable k) - => HashMap k v -> k -> v -(#!) = (HashMap.!) - --- | @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) - --- | 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 - --- | 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]) diff --git a/src/Data/Random/Generics/Internal/Solver.hs b/src/Data/Random/Generics/Internal/Solver.hs deleted file mode 100644 --- a/src/Data/Random/Generics/Internal/Solver.hs +++ /dev/null @@ -1,65 +0,0 @@ --- | Solve systems of equations - -{-# LANGUAGE RecordWildCards #-} -{-# LANGUAGE RankNTypes, FlexibleContexts, TypeFamilies #-} -module Data.Random.Generics.Internal.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) - -> Vector R - -> Maybe (Vector R) -fixedPoint args f = findZero args (liftA2 (V.zipWith (-)) f id) - --- | Assuming @p . f@ is satisfied only for positive values in some interval --- @(0, r]@, find @f r@. -search :: (Double -> a) -> (a -> Bool) -> 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' = 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?" diff --git a/src/Data/Random/Generics/Internal/Types.hs b/src/Data/Random/Generics/Internal/Types.hs deleted file mode 100644 --- a/src/Data/Random/Generics/Internal/Types.hs +++ /dev/null @@ -1,191 +0,0 @@ -{-# LANGUAGE RankNTypes, GADTs, ScopedTypeVariables, ImplicitParams #-} -{-# LANGUAGE TypeOperators, GeneralizedNewtypeDeriving #-} -module Data.Random.Generics.Internal.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 diff --git a/src/Generic/Random/Boltzmann.hs b/src/Generic/Random/Boltzmann.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Boltzmann.hs @@ -0,0 +1,215 @@ +-- | 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, ScopedTypeVariables #-} +{-# LANGUAGE DeriveFunctor, DeriveGeneric, ImplicitParams #-} +{-# LANGUAGE RecordWildCards, DeriveDataTypeable #-} +{-# LANGUAGE TypeFamilies, MultiParamTypeClasses #-} +{-# LANGUAGE TypeApplications #-} +module Generic.Random.Boltzmann 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 Generic.Random.Internal.Common +import Generic.Random.Internal.Solver +import Generic.Random.Internal.Types + +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 @Int)) 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))) + +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]] diff --git a/src/Generic/Random/Data.hs b/src/Generic/Random/Data.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Data.hs @@ -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 Generic.Random.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 Generic.Random.Internal.Data +import Generic.Random.Internal.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 @generic-random@ 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 []) diff --git a/src/Generic/Random/Generic.hs b/src/Generic/Random/Generic.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Generic.hs @@ -0,0 +1,30 @@ +-- | Simple 'GHC.Generics'-based 'arbitrary' generators. +-- +-- Here is an example. Define your type. +-- +-- > data Tree a = Leaf a | Node (Tree a) (Tree a) +-- +-- Derive 'GHC.Generics.Generic'. +-- +-- > deriving 'Generic' -- Turn on the DeriveGeneric extension +-- +-- Pick an arbitrary implementation. +-- +-- > instance Arbitrary a => Arbitrary (Tree a) where +-- > arbitrary = genericArbitraryFrequency [9, 8] +-- +-- @arbitrary :: 'Gen' (Tree a)@ picks a @Leaf@ with probability 9\/17, or a +-- @Node@ with probability 8\/17, and recursively fills their fields with +-- @arbitrary@. + +module Generic.Random.Generic + ( genericArbitrary + , genericArbitraryFrequency + , genericArbitraryFrequency' + , genericArbitrary' + , Nat (..) + , BaseCases' + , BaseCases + ) where + +import Generic.Random.Internal.Generic diff --git a/src/Generic/Random/Internal/Common.hs b/src/Generic/Random/Internal/Common.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Common.hs @@ -0,0 +1,39 @@ +-- | General helper functions + +module Generic.Random.Internal.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 diff --git a/src/Generic/Random/Internal/Data.hs b/src/Generic/Random/Internal/Data.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Data.hs @@ -0,0 +1,146 @@ +{-# LANGUAGE RecordWildCards, DeriveFunctor #-} +module Generic.Random.Internal.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 Generic.Random.Internal.Oracle +import Generic.Random.Internal.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) diff --git a/src/Generic/Random/Internal/Generic.hs b/src/Generic/Random/Internal/Generic.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Generic.hs @@ -0,0 +1,286 @@ +{-# LANGUAGE FlexibleContexts, FlexibleInstances, MultiParamTypeClasses #-} +{-# LANGUAGE TypeApplications, TypeOperators #-} +{-# LANGUAGE DeriveFunctor, GeneralizedNewtypeDeriving #-} +{-# LANGUAGE AllowAmbiguousTypes, ScopedTypeVariables #-} +{-# LANGUAGE DataKinds, KindSignatures #-} +{-# LANGUAGE ConstraintKinds #-} +module Generic.Random.Internal.Generic where + +import Control.Applicative +import Data.Coerce +import GHC.Generics hiding ( S ) +import Test.QuickCheck + +-- * Random generators + +-- | Pick a constructor with uniform probability, and fill its fields +-- recursively. +-- +-- An equivalent definition for @Tree@ is: +-- +-- > genericArbitrary :: Arbitrary a => Gen (Tree a) +-- > genericArbitrary = +-- > oneof +-- > [ Leaf <$> arbitrary -- Uses Arbitrary a +-- > , Node <$> arbitrary <*> arbitrary -- Uses Arbitrary (Tree a) +-- > ] +-- +-- Note that for many types, 'genericArbitrary' tends to produce big values. +-- For instance for @Tree a@ values are finite but the average number of +-- @Leaf@ and @Node@ constructors is infinite. + +genericArbitrary :: (Generic a, GA Unsized (Rep a)) => Gen a +genericArbitrary = ($ repeat 1) . unFreq . fmap to $ ga @Unsized + + +-- | This allows to specify the probability distribution of constructors +-- as a list of weights, in the same order as the data type definition. +-- +-- An equivalent definition for @Tree@ is: +-- +-- > genericArbitraryFrequency :: Arbitrary a => [Int] -> Gen (Tree a) +-- > genericArbitraryFrequency [x, y] = +-- > frequency +-- > [ (x, Leaf <$> arbitrary) +-- > , (y, Node <$> arbitrary <*> arbitrary) +-- > ] + +genericArbitraryFrequency + :: (Generic a, GA Unsized (Rep a)) + => [Int] -- ^ List of weights for every constructor + -> Gen a +genericArbitraryFrequency = unFreq . fmap to $ ga @Unsized + + +-- | The size parameter of 'Gen' is divided among the fields of the chosen +-- constructor. When it reaches zero, the generator selects a finite term +-- whenever it can find any of the given type. +-- +-- The type of 'genericArbitraryFrequency'' has an ambiguous @n@ parameter; it +-- is a type-level natural number of type 'Nat'. That number determines the +-- maximum /depth/ of terms that can be used to end recursion. +-- +-- You'll need the @TypeApplications@ and @DataKinds@ extensions. +-- +-- > genericArbitraryFrequency' @n weights +-- +-- With @n ~ ''Z'@, the generator looks for a simple nullary constructor. If none +-- exist at the current type, as is the case for our @Tree@ type, it carries on +-- as in 'genericArbitraryFrequency'. +-- +-- > genericArbitraryFrequency' @'Z :: Arbitrary a => [Int] -> Gen (Tree a) +-- > genericArbitraryFrequency' @'Z [x, y] = +-- > frequency +-- > [ (x, Leaf <$> arbitrary) +-- > , (y, scale (`div` 2) $ Node <$> arbitrary <*> arbitrary) +-- > ] +-- > -- 2 because Node is 2-ary. +-- +-- Here is another example: +-- +-- > data Tree' = Leaf1 | Leaf2 | Node3 Tree' Tree' Tree' +-- > deriving Generic +-- > +-- > instance Arbitrary Tree' where +-- > arbitrary = genericArbitraryFrequency' @'Z [1, 2, 3] +-- +-- 'genericArbitraryFrequency'' is equivalent to: +-- +-- > genericArbitraryFrequency' @'Z :: [Int] -> Gen Tree' +-- > genericArbitraryFrequency' @'Z [x, y, z] = +-- > sized $ \n -> +-- > if n == 0 then +-- > -- If the size parameter is zero, the non-nullary alternative is discarded. +-- > frequency $ +-- > [ (x, return Leaf1) +-- > , (y, return Leaf2) +-- > ] +-- > else +-- > frequency $ +-- > [ (x, return Leaf1) +-- > , (y, return Leaf2) +-- > , (z, resize (n `div` 3) node) +-- > ] +-- > -- 3 because Node3 is 3-ary +-- > where +-- > node = Node3 <$> arbitrary <*> arbitrary <*> arbitrary +-- +-- To increase the chances of termination when no nullary constructor is directly +-- available, such as in @Tree@, we can pass a larger depth @n@. The effectiveness +-- of this parameter depends on the concrete type the generator is used for. +-- +-- For instance, if we want to generate a value of type @Tree ()@, there is a +-- value of depth 1 (represented by @''S' ''Z'@) that we can use to end +-- recursion: @Leaf ()@. +-- +-- > genericArbitraryFrequency' @('S 'Z) :: [Int] -> Gen (Tree ()) +-- > genericArbitraryFrequency' @('S 'Z) [x, y] = +-- > sized $ \n -> +-- > if n == 0 then +-- > return (Leaf ()) +-- > else +-- > frequency +-- > [ (x, Leaf <$> arbitrary) +-- > , (y, scale (`div` 2) $ Node <$> arbitrary <*> arbitrary) +-- > ] +-- +-- Because the argument of @Tree@ must be inspected in order to discover +-- values of type @Tree ()@, we incur some extra constraints if we want +-- polymorphism. +-- +-- @FlexibleContexts@ and @UndecidableInstances@ are also required. +-- +-- > instance (Arbitrary a, Generic a, BaseCases 'Z (Rep a)) +-- > => Arbitrary (Tree a) where +-- > arbitrary = genericArbitraryFrequency' @('S 'Z) [1, 2] +-- +-- A synonym is provided for brevity. +-- +-- > instance (Arbitrary a, BaseCases' 'Z a) => Arbitrary (Tree a) where +-- > arbitrary = genericArbitraryFrequency' @('S 'Z) [1, 2] + +genericArbitraryFrequency' + :: forall (n :: Nat) a + . (Generic a, GA (Sized n) (Rep a)) + => [Int] -- ^ List of weights for every constructor + -> Gen a +genericArbitraryFrequency' = unFreq . fmap to $ ga @(Sized n) + + +-- | Like 'genericArbitraryFrequency'', but with uniformly distributed +-- constructors. + +genericArbitrary' + :: forall (n :: Nat) a. (Generic a, GA (Sized n) (Rep a)) => Gen a +genericArbitrary' = ($ repeat 1) . unFreq . fmap to $ ga @(Sized n) + + +-- * Internal + +newtype Freq sized a = Freq { unFreq :: [Int] -> Gen a } + deriving Functor + +instance Applicative (Freq sized) where + pure = Freq . pure . pure + Freq f <*> Freq x = Freq (liftA2 (<*>) f x) + +newtype Gen' sized a = Gen' { unGen' :: Gen a } + deriving (Functor, Applicative) + +data Sized :: Nat -> * +data Unsized + +liftGen :: Gen a -> Freq sized a +liftGen = Freq . const + +-- | Generic Arbitrary +class GA sized f where + ga :: Freq sized (f p) + +instance GA sized U1 where + ga = pure U1 + +instance Arbitrary c => GA sized (K1 i c) where + ga = liftGen . fmap K1 $ arbitrary + +instance GA sized f => GA sized (M1 i c f) where + ga = fmap M1 ga + +instance (GASum (Sized n) f, GASum (Sized n) g, BaseCases n f, BaseCases n g) + => GA (Sized n) (f :+: g) where + ga = frequency' gaSum baseCases + where + frequency' :: [Gen' sized a] -> Tagged n [[a]] -> Freq sized a + frequency' as (Tagged a0s) = Freq $ \ws -> + let + units = [(w, elements a0) | (w, a0@(_ : _)) <- zip ws a0s] + in + sized $ \sz -> frequency $ + if sz == 0 && not (null units) then + units + else + [(w, a) | (w, Gen' a) <- zip ws as] + +instance (GASum Unsized f, GASum Unsized g) => GA Unsized (f :+: g) where + ga = frequency' gaSum + where + frequency' :: [Gen' sized a] -> Freq sized a + frequency' as = Freq $ \ws -> frequency + [(w, a) | (w, Gen' a) <- zip ws as] + +instance (GA Unsized f, GA Unsized g) => GA Unsized (f :*: g) where + ga = liftA2 (:*:) ga ga + +instance (GAProduct f, GAProduct g) => GA (Sized n) (f :*: g) where + ga = constScale' a + where + constScale' :: Gen' Unsized a -> Freq (Sized n) a + constScale' = Freq . const . scale (`div` arity) . unGen' + (arity, a) = gaProduct + + +gArbitrarySingle :: forall sized f p . GA sized f => Gen' sized (f p) +gArbitrarySingle = Gen' (unFreq (ga :: Freq sized (f p)) [1]) + + +class GASum sized f where + gaSum :: [Gen' sized (f p)] + +instance (GASum sized f, GASum sized g) => GASum sized (f :+: g) where + gaSum = (fmap . fmap) L1 gaSum ++ (fmap . fmap) R1 gaSum + +instance GA sized f => GASum sized (M1 i c f) where + gaSum = [gArbitrarySingle] + + +class GAProduct f where + gaProduct :: (Int, Gen' Unsized (f p)) + +instance GA Unsized f => GAProduct (M1 i c f) where + gaProduct = (1, gArbitrarySingle) + +instance (GAProduct f, GAProduct g) => GAProduct (f :*: g) where + gaProduct = (m + n, liftA2 (:*:) a b) + where + (m, a) = gaProduct + (n, b) = gaProduct + + +newtype Tagged (a :: Nat) b = Tagged { unTagged :: b } + +-- | Peano-encoded natural numbers. +data Nat = Z | S Nat + +-- | A @BaseCases n ('Rep' a)@ constraint basically provides the list of values +-- of type @a@ with depth at most @n@. +class BaseCases (n :: Nat) f where + baseCases :: Tagged n [[f p]] + +-- | For convenience. +type BaseCases' n a = (Generic a, BaseCases n (Rep a)) + +baseCases' :: forall n f p. BaseCases n f => Tagged n [f p] +baseCases' = (Tagged . concat . unTagged) (baseCases @n) + +instance BaseCases n U1 where + baseCases = Tagged [[U1]] + +instance BaseCases n f => BaseCases n (M1 i c f) where + baseCases = (coerce :: Tagged n [[f p]] -> Tagged n [[M1 i c f p]]) baseCases + +instance BaseCases 'Z (K1 i c) where + baseCases = Tagged [[]] + +instance (Generic c, BaseCases n (Rep c)) => BaseCases ('S n) (K1 i c) where + baseCases = (Tagged . (fmap . fmap) (K1 . to) . unTagged) (baseCases @n) + +instance (BaseCases n f, BaseCases n g) => BaseCases n (f :+: g) where + baseCases = Tagged $ + (fmap . fmap) L1 (unTagged (baseCases @n)) ++ + (fmap . fmap) R1 (unTagged (baseCases @n)) + +instance (BaseCases n f, BaseCases n g) => BaseCases n (f :*: g) where + baseCases = Tagged + [ liftA2 (:*:) + (unTagged (baseCases' @n)) + (unTagged (baseCases' @n)) ] diff --git a/src/Generic/Random/Internal/Oracle.hs b/src/Generic/Random/Internal/Oracle.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Oracle.hs @@ -0,0 +1,499 @@ +{-# LANGUAGE FlexibleContexts, GADTs, RankNTypes, ScopedTypeVariables #-} +{-# LANGUAGE DeriveGeneric, ImplicitParams #-} +{-# LANGUAGE RecordWildCards, DeriveDataTypeable #-} +module Generic.Random.Internal.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 Generic.Random.Internal.Common +import Generic.Random.Internal.Solver +import Generic.Random.Internal.Types + +-- | 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.!) diff --git a/src/Generic/Random/Internal/Solver.hs b/src/Generic/Random/Internal/Solver.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Solver.hs @@ -0,0 +1,66 @@ +-- | Solve systems of equations + +{-# LANGUAGE RecordWildCards #-} +{-# LANGUAGE RankNTypes, FlexibleContexts, TypeFamilies #-} +module Generic.Random.Internal.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?" diff --git a/src/Generic/Random/Internal/Types.hs b/src/Generic/Random/Internal/Types.hs new file mode 100644 --- /dev/null +++ b/src/Generic/Random/Internal/Types.hs @@ -0,0 +1,191 @@ +{-# LANGUAGE RankNTypes, GADTs, ScopedTypeVariables, ImplicitParams #-} +{-# LANGUAGE TypeOperators, GeneralizedNewtypeDeriving #-} +module Generic.Random.Internal.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 diff --git a/test/tree.hs b/test/tree.hs --- a/test/tree.hs +++ b/test/tree.hs @@ -4,7 +4,7 @@ import Data.Foldable import Data.List import Test.QuickCheck -import Data.Random.Generics +import Generic.Random.Data data T = N T T | L deriving (Eq, Ord, Show, Data)