maxent 0.4.0.0 → 0.6.0.0
raw patch · 6 files changed
+107/−298 lines, 6 filesdep +criteriondep +lagrangian
Dependencies added: criterion, lagrangian
Files
- bench/Bench.hs +27/−0
- maxent.cabal +23/−13
- src/Numeric/MaxEnt.hs +12/−9
- src/Numeric/MaxEnt/ConjugateGradient.hs +13/−23
- src/Numeric/MaxEnt/Internal.hs +9/−208
- src/Numeric/MaxEnt/Linear.hs +23/−45
+ bench/Bench.hs view
@@ -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+ ]+ ]
maxent.cabal view
@@ -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.*
src/Numeric/MaxEnt.hs view
@@ -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(..))
src/Numeric/MaxEnt/ConjugateGradient.hs view
@@ -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)
src/Numeric/MaxEnt/Internal.hs view
@@ -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 -
src/Numeric/MaxEnt/Linear.hs view
@@ -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 ------}