diff --git a/bench/Bench.hs b/bench/Bench.hs new file mode 100644 --- /dev/null +++ b/bench/Bench.hs @@ -0,0 +1,27 @@ +{-# LANGUAGE TupleSections, Rank2Types #-} +module Main where +import Numeric.MaxEnt.Internal +import Criterion.Main +import Criterion +import Criterion.Config +import Data.Monoid + +testData :: RealFloat a => ([a], [Constraint a]) +testData = ([1,2,3], [average 1.5]) + +maxent' :: (forall a. RealFloat a => ([a], [Constraint a])) -- ^ A pair of values that the distributions is over and the constraints + -> [Double] +maxent' xs = (\(Right x) -> x) $ maxent 0.00001 xs + + +test1 :: (forall a. RealFloat a => ([a], [Constraint a])) -> Pure +test1 x = nf maxent' x + +myConfig = defaultConfig { cfgReport = Last $ Just "profile.html" , + cfgSamples = Last $ Just 100} + +main = defaultMainWith myConfig (return ()) [ + bgroup "moment" [ + bench "moment1" $ test1 testData + ] + ] diff --git a/maxent.cabal b/maxent.cabal --- a/maxent.cabal +++ b/maxent.cabal @@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.4.0.0 +version: 0.6.0.0 -- A short (one-line) description of the package. synopsis: Compute Maximum Entropy Distributions @@ -33,20 +33,15 @@ . Here is a the example from Probability the Logic of Science . - >>> maxent ([1,2,3], [average 1.5]) - Right [0.61, 0.26, 0.11] + >> maxent 0.00001 [1,2,3] [average 1.5] + > Right [0.61, 0.26, 0.11] . The classic dice example . - >>> maxent ([1,2,3,4,5,6], [average 4.5]) - Right [.05, .07, 0.11, 0.16, 0.23, 0.34] + >> maxent 0.00001 [1,2,3,4,5,6] [average 4.5] + > Right [.05, .07, 0.11, 0.16, 0.23, 0.34] . One can use different constraints besides the average value there. - . - As for why you want to maximize the entropy to find the probability constraint, - I will say this for now. In the case of the average constraint - it is a kin to choosing a integer partition with the most interger compositions. - I doubt that makes any sense, but I will try to explain more with a blog post soon. -- URL for the project homepage or repository. homepage: https://github.com/jfischoff/maxent @@ -80,13 +75,16 @@ 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 -- Other library packages from which modules are imported. build-depends: base ==4.6.*, nonlinear-optimization ==0.3.*, vector ==0.10.*, - ad ==3.4.* + ad ==3.4.*, + lagrangian == 0.5.* -- Directories containing source files. hs-source-dirs: src @@ -101,6 +99,7 @@ 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.*, @@ -109,7 +108,18 @@ 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 @@ -14,12 +14,12 @@ -- -- Here is a the example from Probability the Logic of Science -- --- >>> maxent ([1,2,3], [average 1.5]) +-- >>> maxent 0.00001 [1,2,3] [average 1.5] -- Right [0.61, 0.26, 0.11] -- -- The classic dice example -- --- >>> maxent ([1,2,3,4,5,6], [average 4.5]) +-- >>> maxent 0.00001 [1,2,3,4,5,6] [average 4.5] -- Right [.05, .07, 0.11, 0.16, 0.23, 0.34] -- -- One can use different constraints besides the average value there. @@ -31,28 +31,31 @@ -- module Numeric.MaxEnt ( Constraint, + (.=.), + UU(..), + ExpectationConstraint, ExpectationFunction, - constraint, average, variance, -- ** Classic moment based maxent, -- ** General - GeneralConstraint, - generalMaxent, + general, -- ** Linear LinearConstraints(..), linear, ) where import Numeric.MaxEnt.Internal (Constraint, + (.=.), + UU(..), + ExpectationConstraint, ExpectationFunction, - constraint, average, variance, maxent, - generalMaxent, - GeneralConstraint) -import Numeric.MaxEnt.Linear (linear, LinearConstraints(..)) + general, + 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 @@ -20,31 +20,21 @@ sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c sumWith f xs = sum . zipWith f xs - -toFunction :: (forall a. RealFloat a => [a] -> a) -> Function Simple -toFunction f = VFunction (f . U.toList) - -toGradient :: (forall a. RealFloat a => [a] -> a) -> Gradient Simple -toGradient f = VGradient (U.fromList . grad f . U.toList) - -toDoubleF :: (forall a. RealFloat a => [a] -> a) -> [Double] -> Double -toDoubleF f x = f x - -squaredGrad :: Num a - => (forall s. Mode s => [AD s a] -> AD s a) -> [a] -> a -squaredGrad f vs = sumMap (\x -> x*x) (grad f vs) - -solve :: Double +minimize :: Double -> Int - -> (forall a. RealFloat a => [a] -> a) - -> Either (Result, Statistics) [Double] -solve percision count obj = result where + -> (forall s. Mode s => [AD s Double] -> AD s Double) + -> Either (Result, Statistics) (S.Vector Double) +minimize tolerance count obj = result where guess = U.fromList $ replicate count ((1.0 :: Double) / (fromIntegral count)) - result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) percision guess - (toFunction obj) - (toGradient obj) - Nothing) of - (vs, ToleranceStatisfied, _) -> Right . S.toList $ vs + result = case unsafePerformIO $ + optimize + (defaultParameters { printFinal = False }) + tolerance + guess + (VFunction (lowerFU obj . U.toList)) + (VGradient (U.fromList . grad obj . U.toList)) + Nothing of + (vs, ToleranceStatisfied, _) -> Right vs (_, x, y) -> Left (x, y) 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,210 +1,11 @@ -{-# LANGUAGE TupleSections, Rank2Types #-} -module Numeric.MaxEnt.Internal 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 Numeric.AD.Lagrangian +module Numeric.MaxEnt.Internal ( + module Numeric.MaxEnt.ConjugateGradient, + module Numeric.MaxEnt.General, + module Numeric.MaxEnt.Moment, + module Numeric.MaxEnt.Linear + ) where import Numeric.MaxEnt.ConjugateGradient - -pOfK :: RealFloat a => [a] -> [ExpectationFunction a] -> [a] -> Int -> a -pOfK values fs ls k = exp (negate . sumWith (\l f -> l * f (values !! k)) ls $ fs) / - partitionFunc values fs ls - - - -probs :: (RealFloat b) - => [b] - -> [ExpectationFunction b] - -> [b] - -> [b] -probs values fs ls = map (pOfK values fs ls) [0..length values - 1] - -partitionFunc :: RealFloat a - => [a] - -> [ExpectationFunction a] - -> [a] - -> a -partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | - x <- values, - (f, l) <- zip fs ls] - -objectiveFunc :: RealFloat a - => [a] - -> [ExpectationFunction a] - -> [a] - -> [a] - -> a -objectiveFunc values fs moments ls = log (partitionFunc values fs ls) - + sumWith (*) ls moments - -entropy :: RealFloat a => [a] -> a -entropy = negate . sumMap (\p -> p * log p) - - --- This a functions that takes in a list of values and --- a list of probabilities -type GeneralConstraint a = [a] -> [a] -> a - - -lagrangian :: RealFloat a - => ([a], [GeneralConstraint a], [a]) - -> [a] - -> a -lagrangian (values, fs, constants) lamsAndProbs = result where - result = entropy ps + (sum $ zipWith (*) lams constraints) - constraints = zipWith (-) appliedConstraints constants - appliedConstraints = map (\f -> f values ps) fs - - -- split the args list - ps = take (length values) lamsAndProbs - lams = drop (length values) lamsAndProbs - -generalObjectiveFunc :: RealFloat a => (forall b. RealFloat b => - ([b], [GeneralConstraint b], [b])) - -> [a] - -> a -generalObjectiveFunc params lamsAndProbs = - squaredGrad lang lamsAndProbs where - - lang :: RealFloat a => (forall s. Mode s => [AD s a] -> AD s a) - lang = lagrangian params - - - - --- make a constraint from function and constant -constraint :: RealFloat a => ExpectationFunction a -> a -> Constraint a -constraint = (,) - --- The average constraint -average :: RealFloat a => a -> Constraint a -average m = constraint id m - --- The variance constraint -variance :: RealFloat a => a -> Constraint a -variance sigma = constraint (^(2 :: Int)) sigma - --- | Most general solver --- This is the slowest but most flexible method. Although, I haven't tried using much... - -generalMaxent :: (forall a. RealFloat a => ([a], [(GeneralConstraint a, a)])) -- ^ A pair of values that the distributions is over and the constraints - -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution -generalMaxent params = result where - obj :: RealFloat a => [a] -> a - obj = generalObjectiveFunc objInput - - values :: RealFloat a => [a] - values = fst params - - constraints :: RealFloat a => [(GeneralConstraint a, a)] - constraints = snd params - - fsmoments :: RealFloat a => ([GeneralConstraint a], [a]) - fsmoments = unzip constraints - - --The new constraint is always needed for probability problems - - sumToOne :: RealFloat a => GeneralConstraint a - sumToOne _ xs = sumMap id xs - - gcons' :: RealFloat a => [GeneralConstraint a] - gcons' = fst fsmoments - - gcons :: RealFloat a => [GeneralConstraint a] - gcons = sumToOne : gcons' - - moments :: RealFloat a => [a] - moments = 1 : snd fsmoments - - - objInput :: RealFloat b => ([b], [GeneralConstraint b], [b]) - objInput = (values, gcons, moments) - - fs :: [[Double] -> [Double] -> Double] - fs = fst fsmoments - - --hmm maybe there is a better way to get rid of the defaulting - guess = U.fromList $ replicate - (length fs) (1.0 :: Double) - - result = case unsafePerformIO (optimize defaultParameters 0.00001 guess - (toFunction obj) - (toGradient obj) - Nothing) of - -- Not sure what is supposed to happen here - -- I guess I can plug the lambdas back in - (vs, ToleranceStatisfied, _) -> Right $ (S.toList vs) - (_, x, y) -> Left (x, y) - - --- | Constraint type. A function and the constant it equals. --- --- Think of it as the pair @(f, c)@ in the constraint --- --- @ --- Σ pₐ f(xₐ) = c --- @ --- --- 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 Constraint a = (ExpectationFunction a, a) - --- | A function that takes an index and value and returns a value. --- See 'average' and 'variance' for examples. -type ExpectationFunction a = (a -> a) - --- | Discrete maximum entropy solver where the constraints are all moment constraints. -maxent :: (forall a. RealFloat a => ([a], [Constraint a])) -- ^ A pair of values that the distributions is over and the constraints - -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution -maxent params = result where - obj :: RealFloat a => [a] -> a - obj = uncurry (objectiveFunc values) fsmoments - - values :: RealFloat a => [a] - values = fst params - - constraints :: RealFloat a => [(ExpectationFunction a, a)] - constraints = snd params - - fsmoments :: RealFloat a => ([ExpectationFunction a], [a]) - fsmoments = unzip constraints - - fs :: [Double -> Double] - fs = fst fsmoments - - -- hmm maybe there is a better way to get rid of the defaulting - guess = U.fromList $ replicate - (length fs) (1.0 :: Double) - - {- - result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) - 0.001 guess - (toFunction obj) - (toGradient obj) - Nothing) of - (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs) - (_, x, y) -> Left (x, y) - -} - - result = case unsafePerformIO (optimize (defaultParameters { printFinal = False }) - 0.001 guess - (toFunction obj) - (toGradient obj) - Nothing) of - (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs) - (_, x, y) -> Left (x, y) - - - - - - +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,4 +1,4 @@ -{-# LANGUAGE TupleSections, Rank2Types #-} +{-# LANGUAGE TupleSections, Rank2Types, NoMonomorphismRestriction #-} module Numeric.MaxEnt.Linear where import Numeric.MaxEnt.ConjugateGradient import Numeric.Optimization.Algorithms.HagerZhang05 @@ -11,19 +11,18 @@ import Numeric.AD.Internal.Classes import Data.List (transpose) import Control.Applicative - -probs :: (RealFloat b) => [[b]] -> [b] -> [b] +import qualified Data.Vector.Storable as S + probs matrix ls = map (pOfK matrix ls) [0..length matrix - 1] -pOfK :: RealFloat a => [[a]] -> [a] -> Int -> a pOfK matrix ls k = exp (dot (transpose matrix !! k) ls) / partitionFunc matrix ls -- have the column and rows backwards -partitionFunc :: RealFloat a => [[a]] -> [a] -> a + partitionFunc matrix ws = sum $ [ exp (dot as ws) | as <- transpose matrix] -- This is almost the sam as the objectiveFunc -objectiveFunc :: RealFloat a => [[a]] -> [a]-> [a] -> a + objectiveFunc as moments ls = log (partitionFunc as ls) - dot ls moments data LinearConstraints a = LC { @@ -45,57 +44,36 @@ -- -- >>> 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 :: Double -- ^ Tolerance for the numerical solver - -> (forall a. RealFloat a => LinearConstraints a) -- ^ The matrix A and column vector b - -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution +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 :: RealFloat a => [a] -> a - obj = objectiveFunc (matrix constraints) (output constraints) + obj = objectiveFunc (map (map auto) $ matrix constraints) (map auto $ output constraints) - as :: [[Double]] as = matrix constraints - result = probs as <$> solve tolerance (length as) obj + result = (S.fromList . probs as . S.toList) <$> minimize tolerance count obj ---Just for testing -linear1 :: Double -> (forall a. RealFloat a => LinearConstraints a) -- ^ a matrix A and column vector b - -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution -linear1 precision constraints = result where - obj :: RealFloat a => [a] -> a - obj = objectiveFunc (matrix constraints) (output constraints) - - fs :: [[Double]] - fs = matrix constraints - - guess :: RealFloat a => [a] - guess = replicate - (length fs) ((1.0) / (fromIntegral (length fs))) + count = length $ output constraints + + guess = U.fromList $ replicate + count ((1.0 :: Double) / (fromIntegral count)) + + - result = Right . probs fs . last $ conjugateGradientDescent obj guess -{- - I need to clean this up and submit it - submit the sharpen - - I need to figure out how to make linear2 - what is the - it might help to clean up the notation for the problem - basically I need to symbolically take the derative of the problem - and get rid of the automatic differentation - which I can do hand - - - - --}