diff --git a/bench/Bench.hs b/bench/Bench.hs --- a/bench/Bench.hs +++ b/bench/Bench.hs @@ -5,23 +5,17 @@ import Criterion import Criterion.Config import Data.Monoid +import qualified Data.Vector.Storable as S -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 + 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 @@ -10,7 +10,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 0.6.0.1 +version: 0.6.0.3 -- A short (one-line) description of the package. synopsis: Compute Maximum Entropy Distributions diff --git a/src/Numeric/MaxEnt.hs b/src/Numeric/MaxEnt.hs --- a/src/Numeric/MaxEnt.hs +++ b/src/Numeric/MaxEnt.hs @@ -44,6 +44,8 @@ -- ** Linear LinearConstraints(..), linear, + linear', + linear'' ) where import Numeric.MaxEnt.Internal (Constraint, (.=.), @@ -54,7 +56,9 @@ variance, maxent, general, - linear, + 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 @@ -9,6 +9,7 @@ import Numeric.AD.Types import Numeric.AD.Internal.Classes import Data.List (transpose) +import Control.Arrow (second) dot :: Num a => [a] -> [a] -> a @@ -25,16 +26,25 @@ -> (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)) + guess = U.fromList $ 1 : replicate (count - 1) 0 result = case unsafePerformIO $ optimize - (defaultParameters { printFinal = False }) + (defaultParameters { printFinal = False, + --printParams = True, + --maxit = 10000, + --verbosity = VeryVerbose, + initialStep = Just 0.1 + --lineSearch = ApproximateWolfe, + --debugTol = Just 0.01, + --nanRho = 1.3, + --estimateError = RelativeEpsilon 0.1 + --lbfgs = False + }) tolerance guess (VFunction (lowerFU obj . U.toList)) (VGradient (U.fromList . grad obj . U.toList)) - Nothing of + (Just $ VCombined (second U.fromList . grad' obj . U.toList)) of (vs, ToleranceStatisfied, _) -> Right vs (_, x, y) -> Left (x, y) 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,6 +1,6 @@ {-# LANGUAGE TupleSections, Rank2Types, NoMonomorphismRestriction #-} module Numeric.MaxEnt.Linear where -import Numeric.MaxEnt.ConjugateGradient +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 @@ -13,17 +13,17 @@ import Control.Applicative import qualified Data.Vector.Storable as S -probs matrix ls = map (pOfK matrix ls) [0..length matrix - 1] - -pOfK matrix ls k = exp (dot (transpose matrix !! k) ls) / partitionFunc matrix ls - --- have the column and rows backwards +multMV mat vec = map (\row -> dot row vec) mat + +probs matrix ls = result where + norm = partitionFunc matrix ls + result = map (\x -> exp x / norm ) $ (transpose matrix) `multMV` ls -partitionFunc matrix ws = sum $ [ exp (dot as ws) | as <- transpose matrix] +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 as moments ls = (log (partitionFunc as ls) - dot ls moments) data LinearConstraints a = LC { matrix :: [[a]], @@ -58,26 +58,46 @@ linear tolerance constraints = result where obj = objectiveFunc (map (map auto) $ matrix constraints) (map auto $ output constraints) - as = matrix constraints - - result = (S.fromList . probs as . S.toList) <$> minimize tolerance count obj - + as = matrix constraints count = length $ output constraints - guess = U.fromList $ replicate - count ((1.0 :: Double) / (fromIntegral count)) - + 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) + as = matrix constraints + count = length $ output constraints + guess = 1 : replicate (count - 1) 0 + result = map (probs as) . gradientDescent obj $ guess + + +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) + 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] 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 @@ -54,22 +54,21 @@ variance :: Num a => a -> ExpectationConstraint a variance sigma = (^(2 :: Int)) .=. sigma -pOfK :: [Double] -> [ExpectationFunction Double] -> S.Vector Double -> Int -> Double -pOfK values fs ls k = - exp (negate . sum . zipWith (\l f -> l * f (values !! k)) lsList $ fs) / - (partitionFunc values fs lsList) where - lsList = S.toList ls +--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 = S.map (pOfK values fs ls) . S.enumFromN 0 $ length 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 -partitionFunc values fs ls = sum $ [ exp ((-l) * f x) | - x <- values, - (f, l) <- zip fs ls] +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 (\x y -> x * y) ls moments) + log (partitionFunc values fs ls) + (sum $ zipWith (*) ls moments) --- My thoughts are that I should maybe split this up -- | Discrete maximum entropy solver where the constraints are all moment constraints. maxent :: Double -- ^ Tolerance for the numerical solver diff --git a/tests/Main.hs b/tests/Main.hs --- a/tests/Main.hs +++ b/tests/Main.hs @@ -4,6 +4,18 @@ import Test.Framework.Providers.QuickCheck2 (testProperty) import Numeric.MaxEnt import LinearTesting +import qualified Data.Vector.Storable as S +import Data.List (intersperse) +import Numeric + +un (Right x) = S.toList x + +pg (Right x) = putStrLn . concat . intersperse ", " . + map (($ "") . showFFloat (Just 3)) . S.toList . fst $ x +pg (Left x) = print x + +pm (Right x) = putStrLn . concat . intersperse ", " . map (($ "") . showFFloat (Just 3)) . S.toList $ x +pm (Left x) = print x -- One test is that a list convolved with a -- guassian and then infered with a gaussian should have less entropy then it