diff --git a/bench/Bench.hs b/bench/Bench.hs deleted file mode 100644 --- a/bench/Bench.hs +++ /dev/null @@ -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]) - ] - ] diff --git a/maxent.cabal b/maxent.cabal --- a/maxent.cabal +++ b/maxent.cabal @@ -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.* - - - diff --git a/src/Numeric/MaxEnt.hs b/src/Numeric/MaxEnt.hs --- a/src/Numeric/MaxEnt.hs +++ b/src/Numeric/MaxEnt.hs @@ -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(..)) - - - - - diff --git a/src/Numeric/MaxEnt/ConjugateGradient.hs b/src/Numeric/MaxEnt/ConjugateGradient.hs --- a/src/Numeric/MaxEnt/ConjugateGradient.hs +++ b/src/Numeric/MaxEnt/ConjugateGradient.hs @@ -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 diff --git a/src/Numeric/MaxEnt/General.hs b/src/Numeric/MaxEnt/General.hs --- a/src/Numeric/MaxEnt/General.hs +++ b/src/Numeric/MaxEnt/General.hs @@ -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 diff --git a/src/Numeric/MaxEnt/Internal.hs b/src/Numeric/MaxEnt/Internal.hs --- a/src/Numeric/MaxEnt/Internal.hs +++ b/src/Numeric/MaxEnt/Internal.hs @@ -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 diff --git a/src/Numeric/MaxEnt/Linear.hs b/src/Numeric/MaxEnt/Linear.hs --- a/src/Numeric/MaxEnt/Linear.hs +++ b/src/Numeric/MaxEnt/Linear.hs @@ -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] +-- ) diff --git a/src/Numeric/MaxEnt/Moment.hs b/src/Numeric/MaxEnt/Moment.hs --- a/src/Numeric/MaxEnt/Moment.hs +++ b/src/Numeric/MaxEnt/Moment.hs @@ -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 diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -1,4 +1,5 @@ module Main where + import Test.Framework (defaultMain, testGroup, defaultMainWithArgs) import Test.Framework.Providers.HUnit import Test.Framework.Providers.QuickCheck2 (testProperty)