module 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)
sumWith :: Num c => (a -> b -> c) -> [a] -> [b] -> c
sumWith f xs = sum . zipWith f xs
pOfK :: Floating a => [a] -> [ExpectationFunction a] -> [a] -> Int -> a
pOfK values fs ls k = exp (negate . sumWith (\l f -> l * f k (values !! k)) ls $ fs) /
partitionFunc values fs ls
probs :: Floating b
=> [b]
-> [ExpectationFunction b]
-> [b]
-> [b]
probs values fs ls = map (pOfK values fs ls) [0..length values 1]
partitionFunc :: Floating a
=> [a]
-> [ExpectationFunction a]
-> [a]
-> a
partitionFunc values fs ls = sum $ [ exp ((l) * f i x) |
(i, x) <- zip [0..] values,
(f, l) <- zip fs ls]
objectiveFunc :: Floating a
=> [a]
-> [ExpectationFunction a]
-> [a]
-> [a]
-> a
objectiveFunc values fs moments ls = log (partitionFunc values fs ls)
+ sumWith (*) ls moments
toFunction :: (forall a. Floating a => [a] -> a) -> Function Simple
toFunction f = VFunction (f . U.toList)
toGradient :: (forall a. Floating a => [a] -> a) -> Gradient Simple
toGradient f = VGradient (U.fromList . grad f . U.toList)
toDoubleF :: (forall a. Floating a => [a] -> a) -> [Double] -> Double
toDoubleF f x = f x
type Constraint a = (ExpectationFunction a, a)
type ExpectationFunction a = (Int -> a -> a)
constraint :: Floating a => ExpectationFunction a -> a -> Constraint a
constraint = (,)
average :: Floating a => a -> Constraint a
average m = constraint (const id) m
variance :: Floating a => a -> Constraint a
variance sigma = constraint (const (^(2 :: Int))) sigma
maxent :: (forall a. Floating a => ([a], [Constraint a]))
-> Either (Result, Statistics) [Double]
maxent params = result where
obj :: Floating a => [a] -> a
obj = uncurry (objectiveFunc values) fsmoments
values :: Floating a => [a]
values = fst params
constraints :: Floating a => [(ExpectationFunction a, a)]
constraints = snd params
fsmoments :: Floating a => ([ExpectationFunction a], [a])
fsmoments = unzip constraints
fs :: [Int -> Double -> Double]
fs = fst fsmoments
guess = U.fromList $ replicate
(length fs) (1.0 :: Double)
result = case unsafePerformIO (optimize defaultParameters 0.00001 guess
(toFunction obj)
(toGradient obj)
Nothing) of
(vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs)
(_, x, y) -> Left (x, y)