maxent 0.6.0.4 → 0.7
raw patch · 9 files changed
+164/−236 lines, 9 filesdep −criteriondep ~QuickCheckdep ~addep ~base
Dependencies removed: criterion
Dependency ranges changed: QuickCheck, ad, base, hmatrix, lagrangian
Files
- bench/Bench.hs +0/−21
- maxent.cabal +19/−39
- src/Numeric/MaxEnt.hs +6/−10
- src/Numeric/MaxEnt/ConjugateGradient.hs +25/−13
- src/Numeric/MaxEnt/General.hs +8/−14
- src/Numeric/MaxEnt/Internal.hs +3/−2
- src/Numeric/MaxEnt/Linear.hs +73/−80
- src/Numeric/MaxEnt/Moment.hs +29/−57
- tests/Main.hs +1/−0
− bench/Bench.hs
@@ -1,21 +0,0 @@-{-# LANGUAGE TupleSections, Rank2Types #-}-module Main where-import Numeric.MaxEnt.Internal-import Criterion.Main-import Criterion-import Criterion.Config-import Data.Monoid-import qualified Data.Vector.Storable as S- ---myConfig = defaultConfig { cfgReport = Last $ Just "profile.html" ,- cfgSamples = Last $ Just 100}--main = defaultMainWith myConfig (return ()) [- bgroup "linear" [- bench "linear1" $ nf ((\(Right x) -> x) . linear 3.0e-17) (LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298]),- bench "linear'" $ nf ((\(Right x) -> x) . linear 3.0e-17) (LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298]),- bench "linear''" $ nf ((\(Right x) -> x) . linear 3.0e-17) (LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298])- ]- ]
maxent.cabal view
@@ -4,13 +4,7 @@ -- The name of the package. name: maxent --- The package version. See the Haskell package versioning policy (PVP) --- for standards guiding when and how versions should be incremented.--- http://www.haskell.org/haskellwiki/Package_versioning_policy--- PVP summary: +-+------- breaking API changes--- | | +----- non-breaking API additions--- | | | +--- code changes with no API change-version: 0.6.0.4+version: 0.7 -- A short (one-line) description of the package. synopsis: Compute Maximum Entropy Distributions@@ -53,7 +47,7 @@ license-file: LICENSE -- The package author(s).-author: Jonathan Fischoff+author: (c) Jonathan Fischoff 2012-2014, (c) Eric Pashman 2014 -- An email address to which users can send suggestions, bug reports, and -- patches.@@ -71,22 +65,23 @@ library+ ghc-options: -Wall -- Modules exported by the library. exposed-modules: Numeric.MaxEnt -- Modules included in this library but not exported.- other-modules: Numeric.MaxEnt.Internal, - Numeric.MaxEnt.Linear, - Numeric.MaxEnt.ConjugateGradient,+ other-modules: Numeric.MaxEnt.Internal,+ Numeric.MaxEnt.Linear,+ Numeric.MaxEnt.ConjugateGradient Numeric.MaxEnt.Moment, Numeric.MaxEnt.General -- Other library packages from which modules are imported.- build-depends: base ==4.6.*,- nonlinear-optimization ==0.3.*,- vector ==0.10.*, - ad ==3.4.*,- lagrangian == 0.5.*+ build-depends: base >=4.5 && < 5,+ nonlinear-optimization == 0.3.*,+ vector == 0.10.*, + ad >= 4 && < 5,+ lagrangian == 0.6.* -- Directories containing source files. hs-source-dirs: src@@ -96,32 +91,17 @@ hs-source-dirs: src, tests type: exitcode-stdio-1.0 main-is: Main.hs- build-depends: base ==4.6.*,+ build-depends: base >=4.5 && < 5, nonlinear-optimization ==0.3.*, vector ==0.10.*, - ad ==3.4.*,- hmatrix ==0.14.*,- lagrangian == 0.5.*,- QuickCheck == 2.5.*,- test-framework-quickcheck2 ==0.3.*,- test-framework-quickcheck2 ==0.3.*,- test-framework-hunit ==0.3.*,+ ad >= 4 && < 5,+ hmatrix >= 0.14 && < 0.17,+ lagrangian == 0.6.*,+ QuickCheck,+ test-framework-quickcheck2 == 0.3.*,+ test-framework-quickcheck2 == 0.3.*,+ test-framework-hunit == 0.3.*, test-framework == 0.8.* default-language: Haskell2010 -Benchmark bench- default-language: Haskell2010- hs-source-dirs: src, bench- type: exitcode-stdio-1.0- main-is: Bench.hs- build-depends: base ==4.6.*,- nonlinear-optimization ==0.3.*,- vector ==0.10.*, - ad ==3.4.*,- hmatrix ==0.14.*,- lagrangian == 0.5.*,- criterion == 0.6.*---
src/Numeric/MaxEnt.hs view
@@ -32,9 +32,7 @@ module Numeric.MaxEnt ( Constraint, (.=.),- UU(..), ExpectationConstraint,- ExpectationFunction, average, variance, -- ** Classic moment based@@ -43,21 +41,19 @@ general, -- ** Linear LinearConstraints(..),- linear+ linear,+ linear',+ linear'' ) where+ import Numeric.MaxEnt.Internal (Constraint, (.=.),- UU(..), ExpectationConstraint,- ExpectationFunction, average, variance, maxent, general, linear,+ linear',+ linear'', LinearConstraints(..))-----
src/Numeric/MaxEnt/ConjugateGradient.hs view
@@ -1,29 +1,41 @@-{-# LANGUAGE TupleSections, Rank2Types #-}+{-# LANGUAGE Rank2Types #-}++--------------------------------------------------------------------------------+-- This module is updated to work with version 4.* of `Numeric.AD`, but it is+-- now provides funcationality only to `Numeric.MaxEnt.Linear`. Formerly,+-- `Numeric.MaxEnt.Moment` used the `minimize` function defined here, but I+-- rewrote that module to use the `general` function defined in+-- `Numeric.MaxEnt.General`, which in turn uses `maximize` from the+-- `Numeric.AD.Lagrangian`.+--+-- I intend to rewrite `Numeric.MaxEnt.Linear` so that it no longer relies on+-- this module either, with would leave it unused. -- E.P.+--------------------------------------------------------------------------------+ module Numeric.MaxEnt.ConjugateGradient where-import Numeric.Optimization.Algorithms.HagerZhang05++import Control.Arrow (second)+ import qualified Data.Vector.Unboxed as U import qualified Data.Vector.Storable as S-import Numeric.AD-import GHC.IO (unsafePerformIO)-import Data.Traversable-import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Data.List (transpose)-import Control.Arrow (second) +import GHC.IO (unsafePerformIO) +import Numeric.Optimization.Algorithms.HagerZhang05+import Numeric.AD+ dot :: Num a => [a] -> [a] -> a-dot x y = sum . zipWith (*) x $ y+dot xs ys = sum $ zipWith (*) xs ys sumMap :: Num b => (a -> b) -> [a] -> b sumMap f = sum . map f sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c -sumWith f xs = sum . zipWith f xs+sumWith f xs ys = sum $ zipWith f xs ys minimize :: Double -> Int- -> (forall s. Mode s => [AD s Double] -> AD s Double) + -> (forall a. (Floating a) => [a] -> a) -> Either (Result, Statistics) (S.Vector Double) minimize tolerance count obj = result where guess = U.fromList $ 1 : replicate (count - 1) 0@@ -43,7 +55,7 @@ }) tolerance guess - (VFunction (lowerFU obj . U.toList))+ (VFunction (obj . U.toList)) (VGradient (U.fromList . grad obj . U.toList)) (Just $ VCombined (second U.fromList . grad' obj . U.toList)) of (vs, ToleranceStatisfied, _) -> Right vs
src/Numeric/MaxEnt/General.hs view
@@ -1,19 +1,15 @@ {-# LANGUAGE TupleSections, Rank2Types #-}+ module Numeric.MaxEnt.General ( Constraint, general ) where-import Numeric.Optimization.Algorithms.HagerZhang05-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Storable as S-import Numeric.AD-import GHC.IO (unsafePerformIO)-import Data.Traversable-import Numeric.AD.Types-import Numeric.AD.Internal.Tower-import Numeric.AD.Internal.Classes-import Data.List (transpose)+ import Control.Applicative++import qualified Data.Vector.Storable as S++import Numeric.Optimization.Algorithms.HagerZhang05 (Result, Statistics) import Numeric.AD.Lagrangian entropy :: Floating a => [a] -> a@@ -25,11 +21,9 @@ -- ^ Tolerance for the numerical solver -> Int -- ^ the count of probabilities- -> [Constraint Double]+ -> [Constraint] -- ^ constraints -> Either (Result, Statistics) (S.Vector Double) -- ^ Either the a discription of what wrong or the probability distribution general tolerance count constraints = - fst <$> maximize tolerance entropy ((sum <=> 1.0) : constraints) count- - + fst <$> maximize entropy ((sum <=> 1) : constraints) tolerance count
src/Numeric/MaxEnt/Internal.hs view
@@ -1,10 +1,11 @@ module Numeric.MaxEnt.Internal (- module Numeric.MaxEnt.ConjugateGradient,+ --module Numeric.MaxEnt.ConjugateGradient, module Numeric.MaxEnt.General, module Numeric.MaxEnt.Moment, module Numeric.MaxEnt.Linear ) where-import Numeric.MaxEnt.ConjugateGradient++--import Numeric.MaxEnt.ConjugateGradient import Numeric.MaxEnt.General import Numeric.MaxEnt.Moment import Numeric.MaxEnt.Linear
src/Numeric/MaxEnt/Linear.hs view
@@ -1,36 +1,41 @@-{-# LANGUAGE TupleSections, Rank2Types, NoMonomorphismRestriction #-}+{-# LANGUAGE FlexibleContexts, Rank2Types, NoMonomorphismRestriction,+ StandaloneDeriving #-}+ module Numeric.MaxEnt.Linear where-import Numeric.MaxEnt.ConjugateGradient (minimize, dot)-import Numeric.Optimization.Algorithms.HagerZhang05-import qualified Data.Vector.Unboxed as U-import qualified Data.Vector.Storable as S-import Numeric.AD-import GHC.IO (unsafePerformIO)-import Data.Traversable-import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Data.List (transpose)+ import Control.Applicative++import Data.List (transpose) import qualified Data.Vector.Storable as S++import Numeric.MaxEnt.ConjugateGradient (minimize, dot)+import Numeric.Optimization.Algorithms.HagerZhang05 (Result, Statistics)+import Numeric.AD +multMV :: (Num a) => [[a]] -> [a] -> [a] multMV mat vec = map (\row -> dot row vec) mat +probs :: (Floating a) => [[a]] -> [a] -> [a] probs matrix ls = result where norm = partitionFunc matrix ls result = map (\x -> exp x / norm ) $ (transpose matrix) `multMV` ls +partitionFunc :: (Floating a) => [[a]] -> [a] -> a partitionFunc matrix ws = sum . map exp . multMV (transpose matrix) $ ws -- This is almost the sam as the objectiveFunc -objectiveFunc as moments ls = (log (partitionFunc as ls) - dot ls moments)+objectiveFunc :: (Floating a) => [[a]] -> [a] -> [a] -> a+objectiveFunc as moments ls = log $ partitionFunc as ls - dot ls moments -data LinearConstraints a = LC {- matrix :: [[a]], - output :: [a]- }- deriving (Show, Eq)+data LinearConstraints = LC+ { unLC :: forall a. (Floating a) => ([[a]], [a]) } +-- These instances default the underlying numeric type of `LC` to `Double`,+-- which may be problematic for some usages.+deriving instance Eq LinearConstraints+deriving instance Show LinearConstraints+ -- | This is for the linear case Ax = b -- @x@ in this situation is the vector of probablities. -- @@ -42,72 +47,60 @@ -- -- Now if we were given just the convolution and the output, we can use 'linear' to infer the input. -- --- >>> linear 3.0e-17 $ LC [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] [0.276, 0.426, 0.298]--- Right [0.20000000000000004,0.4999999999999999,0.3]--- --- I fell compelled to point out that we could also just invert the original convolution --- matrix. Supposedly using maxent can reduce errors from noise if the convolution --- matrix is not properly estimated.--- +-- >>> linear 3.0e-17 $ LC ([[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]], [0.276, 0.426, 0.298])+-- Right (fromList [0.2000000000000001,0.49999999999999983,0.30000000000000004])+--+-- I fell compelled to point out that we could also just invert the original+-- convolution matrix. Supposedly using maxent can reduce errors from noise if+-- the convolution matrix is not properly estimated. linear :: Double - -- ^ Tolerance for the numerical solver- -> LinearConstraints Double- -- ^ The matrix A and column vector b- -> Either (Result, Statistics) (S.Vector Double)- -- ^ Either the a discription of what wrong or the probability distribution -linear tolerance constraints = result where- obj = objectiveFunc (map (map auto) $ matrix constraints) (map auto $ output constraints) -- as = matrix constraints- count = length $ output constraints - - result = (S.fromList . probs as . S.toList) <$> minimize tolerance count obj- --linear' :: LinearConstraints Double- -- ^ The matrix A and column vector b- -> [[Double]]- -- ^ Either the a discription of what wrong or the probability distribution-linear' constraints = result where- obj = objectiveFunc (map (map auto) $ matrix constraints) (map auto $ output constraints) + -- ^ Tolerance for the numerical solver+ -> LinearConstraints+ -- ^ The matrix A and column vector b+ -> Either (Result, Statistics) (S.Vector Double)+ -- ^ Either a description of what went wrong or the probability+ -- distribution +linear tolerance constraints =+ let (matrix, output) = unLC constraints+ obj = objectiveFunc matrix output + n = length output+ in (S.fromList . probs matrix . S.toList) <$> minimize tolerance n obj - as = matrix constraints- count = length $ output constraints - guess = 1 : replicate (count - 1) 0+--------------------------------------------------------------------------------+-- I updated everything below to work with the new types, but it's not clear to +-- me what it's for. -- EP+-------------------------------------------------------------------------------- - result = map (probs as) . gradientDescent obj $ guess- +linear' :: (Floating a, Ord a)+ => LinearConstraints+ -- ^ The matrix A and column vector b+ -> [[a]]+ -- ^ Either a description of what went wrong or the probability+ -- distribution+linear' constraints =+ let (matrix, output) = unLC constraints+ obj = objectiveFunc matrix output+ guess = 1 : replicate (length output - 1) 0+ in map (probs matrix) . gradientDescent obj $ guess -linear'' :: LinearConstraints Double+linear'' :: (Floating a, Ord a)+ => LinearConstraints -- ^ The matrix A and column vector b- -> [[Double]]- -- ^ Either the a discription of what wrong or the probability distribution-linear'' constraints = result where- obj = objectiveFunc (map (map auto) $ matrix constraints) (map auto $ output constraints) -- as = matrix constraints- count = length $ output constraints - guess = 1 : replicate (count - 1) 0-- result = map (probs as) . conjugateGradientDescent obj $ guess--test1 = LC - [[0.892532,0.003851,0.063870,0.001593,0.038155],- [0.237713,0.111149,0.326964,0.271535,0.052639],- [0.133708,0.788233,0.051543,0.003976,0.022539],- [0.238064,0.263171,0.112279,0.270452,0.116034],- [0.844155,0.011312,0.001470,0.001826,0.141237]]- [0.246323,0.235600,0.071699,0.211339,0.238439]------ - ----+ -> [[a]]+ -- ^ Either a description of what went wrong or the probability+ -- distribution+linear'' constraints =+ let (matrix, output) = unLC constraints+ obj = objectiveFunc matrix output + guess = 1 : replicate (length output - 1) 0+ in map (probs matrix) . conjugateGradientDescent obj $ guess - +--test1 = LC ( [ [0.892532,0.003851,0.063870,0.001593,0.038155]+-- , [0.237713,0.111149,0.326964,0.271535,0.052639]+-- , [0.133708,0.788233,0.051543,0.003976,0.022539]+-- , [0.238064,0.263171,0.112279,0.270452,0.116034]+-- , [0.844155,0.011312,0.001470,0.001826,0.141237]+-- ]+-- ,+-- [0.246323,0.235600,0.071699,0.211339,0.238439]+-- )
src/Numeric/MaxEnt/Moment.hs view
@@ -1,27 +1,19 @@-{-# LANGUAGE TupleSections, Rank2Types, NoMonomorphismRestriction #-}+{-# LANGUAGE Rank2Types, NoMonomorphismRestriction #-}+ module Numeric.MaxEnt.Moment ( ExpectationConstraint, (.=.),- ExpectationFunction, average, variance,- maxent,- UU(..)+ maxent ) where-import Numeric.Optimization.Algorithms.HagerZhang05-import qualified Data.Vector.Unboxed as U+ import qualified Data.Vector.Storable as S-import Numeric.AD-import GHC.IO (unsafePerformIO)-import Data.Traversable-import Numeric.AD.Types-import Numeric.AD.Internal.Classes-import Data.List (transpose)-import Control.Applicative-import Numeric.MaxEnt.ConjugateGradient-import Data.List (foldl')---import Data.Vector +import Numeric.Optimization.Algorithms.HagerZhang05 (Result, Statistics)+import Numeric.AD.Lagrangian+import Numeric.MaxEnt.General+ -- | Constraint type. A function and the constant it equals. -- -- Think of it as the pair @(f, c)@ in the constraint @@ -33,60 +25,40 @@ -- such that we are summing over all values . -- -- For example, for a variance constraint the @f@ would be @(\\x -> x*x)@ and @c@ would be the variance.-type ExpectationConstraint a = (UU a, a)+newtype ExpectationConstraint = ExpCon+ { unExpCon :: forall a. (Floating a) => [a] -> ([a] -> a, a) } ----infixr 1 .=.-(.=.) :: (forall s. Mode s => AD s a -> AD s a) -> a -> ExpectationConstraint a-f .=. c = (UU f, c) --- | A function that takes an index and value and returns a value.--- See 'average' and 'variance' for examples.-type ExpectationFunction a = (a -> a)+infixr 1 .=.+(.=.) :: (forall a. (Floating a) => a -> a)+ -> (forall b. (Floating b) => b)+ -> ExpectationConstraint+f .=. c = ExpCon $ \vals -> (sum .zipWith (*) vals . map f , c) -newtype UU a = UU {unUU :: forall s. Mode s => ExpectationFunction (AD s a) }+expCon2Con :: (forall a. (Floating a) => [a])+ -> ExpectationConstraint+ -> Constraint+expCon2Con vals expCon = f <=> c where+ (f, c) = unExpCon expCon vals -- The average constraint-average :: Num a => a -> ExpectationConstraint a+average :: (forall a. (Floating a) => a) -> ExpectationConstraint average m = id .=. m -- The variance constraint-variance :: Num a => a -> ExpectationConstraint a+variance :: (forall a. (Floating a) => a) -> ExpectationConstraint variance sigma = (^(2 :: Int)) .=. sigma ---partialPart' ls fs x = exp . negate . S.sum . S.zipWith (\l f -> l * f x) ls $ fs---partitionFunc' values fs ls = S.sum . S.map (partialPart' ls fs) $ values--probs values fs ls = result where- lsList = S.toList ls- norm = partitionFunc values fs lsList- result = S.map (\x -> partialPart lsList fs x / norm) $ S.fromList values --partialPart ls fs x = exp . negate . sum . zipWith (\l f -> l * f x) ls $ fs--partitionFunc values fs ls = sum . map (partialPart ls fs) $ values--objectiveFunc fs moments values ls = - log (partitionFunc values fs ls) + (sum $ zipWith (*) ls moments)---- | Discrete maximum entropy solver where the constraints are all moment constraints. +-- | Discrete maximum entropy solver where the constraints are all moment+-- constraints. maxent :: Double -- ^ Tolerance for the numerical solver- -> [Double]+ -> (forall a. (Floating a) => [a]) -- ^ values that the distributions is over- -> [ExpectationConstraint Double]+ -> [ExpectationConstraint] -- ^ The constraints -> Either (Result, Statistics) (S.Vector Double) -- ^ Either the a discription of what wrong or the probability distribution -maxent tolerance values constraints = result where- obj = objectiveFunc (map unUU fs') (map auto moments) (map auto values)- - count = length fs- - (fs', moments) = unzip constraints - - fs = map (\x -> lowerUU $ unUU x) fs'- - guess = U.fromList $ replicate count (1.0 / fromIntegral count :: Double) - - result = probs values fs <$> minimize tolerance count obj+maxent tolerance values expConstraints = general tolerance n constraints where+ constraints = map (expCon2Con values) expConstraints + n = length values
tests/Main.hs view
@@ -1,4 +1,5 @@ module Main where+ import Test.Framework (defaultMain, testGroup, defaultMainWithArgs) import Test.Framework.Providers.HUnit import Test.Framework.Providers.QuickCheck2 (testProperty)