packages feed

maxent 0.3.1.1 → 0.4.0.0

raw patch · 8 files changed

+482/−363 lines, 8 filesdep +QuickCheckdep +hmatrixdep +test-frameworkdep ~ad

Dependencies added: QuickCheck, hmatrix, test-framework, test-framework-hunit, test-framework-quickcheck2

Dependency ranges changed: ad

Files

maxent.cabal view
@@ -10,7 +10,7 @@ -- PVP summary:      +-+------- breaking API changes --                   | | +----- non-breaking API additions --                   | | | +--- code changes with no API change-version:             0.3.1.1+version:             0.4.0.0  -- A short (one-line) description of the package. synopsis:            Compute Maximum Entropy Distributions@@ -33,18 +33,12 @@  .  Here is a the example from Probability the Logic of Science  .- @-  maxent ([1,2,3], [average 1.5])- @- .+ >>> maxent ([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 ([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.  @@ -78,22 +72,44 @@ build-type:          Simple  -- Constraint on the version of Cabal needed to build this package.-cabal-version:       >=1.8+cabal-version:       >=1.12   library   -- Modules exported by the library.-  exposed-modules:     MaxEnt+  exposed-modules:     Numeric.MaxEnt      -- Modules included in this library but not exported.-  other-modules:     MaxEnt.Internal+  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.2.*+                 ad ==3.4.*      -- Directories containing source files.   hs-source-dirs:      src+  default-language: Haskell2010+  +Test-Suite tests+  hs-source-dirs:      src, tests+  type:       exitcode-stdio-1.0+  main-is:    Main.hs+  build-depends: base ==4.6.*,+                 nonlinear-optimization ==0.3.*,+                 vector ==0.10.*, +                 ad ==3.4.*,+                 hmatrix ==0.14.*,+                 QuickCheck == 2.5.*,+                 test-framework-quickcheck2 ==0.3.*,+                 test-framework-quickcheck2 ==0.3.*,+                 test-framework-hunit ==0.3.*,+                 test-framework == 0.8.*+                    +  default-language: Haskell2010++++
− src/MaxEnt.hs 
@@ -1,55 +0,0 @@--- |--- The maximum entropy method, or MAXENT, is variational approach for computing probability --- distributions given a list of moment, or expected value, constraints.--- --- Here are some links for background info.--- A good overview of applications:--- <http://cmm.cit.nih.gov/maxent/letsgo.html>--- On the idea of maximum entropy in general: --- <http://en.wikipedia.org/wiki/Principle_of_maximum_entropy>---  --- --- Use this package to compute discrete maximum entropy distributions over a list of values and--- list of constraints.--- --- 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]--- --- 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]--- --- 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.--- -module MaxEnt (-    Constraint,-    ExpectationFunction,-    constraint,-    average,-    variance,-    maxent,-    maxentLinear,-    generalMaxent-) where-import MaxEnt.Internal (Constraint,-                        ExpectationFunction,-                        constraint,-                        average,-                        variance,-                        maxent,-                        maxentLinear,-                        generalMaxent)
− src/MaxEnt/Internal.hs 
@@ -1,295 +0,0 @@-{-# LANGUAGE TupleSections, Rank2Types #-}-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)-import Data.Traversable-import Numeric.AD.Types-import Numeric.AD.Internal.Classes---import Numeric.AD.Lagrangian--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--pOfK :: Floating 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 --pOfKLinear :: Floating a => [[a]] -> [a] -> Int -> a-pOfKLinear matrix ls k = -    exp (dot (matrix !! k) ls) / -        partitionFuncLinear matrix 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 x) | -                                x <- 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---dot :: Num a => [a] -> [a] -> a-dot x y = sum . zipWith (*) x $ y--partitionFuncLinear :: Floating a-             => [[a]]-             -> [a] -             -> a-partitionFuncLinear matrix ws = sum $ [ exp (dot as ws) | -                               as <- matrix]---- This is almost the sam as the objectiveFunc                                   -objectiveFuncLinear :: Floating a-             => [a] -             -> [[a]] -             -> [a] -             -> [a] -             -> a-objectiveFuncLinear values as moments ls = -    log (partitionFuncLinear as ls) - dot ls moments---linProbs :: (Floating b)-      => [[b]] -      -> [b] -      -> [b]    -linProbs matrix ls = -    map (pOfKLinear matrix ls) [0..length matrix - 1]---entropy :: Floating 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 :: Floating 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--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)--generalObjectiveFunc :: Floating a => (forall b. Floating b => -             ([b], [GeneralConstraint b], [b]))-             -> [a] -             -> a-generalObjectiveFunc params lamsAndProbs = -    squaredGrad lang lamsAndProbs  where-        -    lang :: Floating a => (forall s. Mode s => [AD s a] -> AD s a)-    lang = lagrangian params---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 ---- | Constraint type. A function and the constant it equals.--- ---   Think of it as the pair @(f, c)@ in the constraint ------ @---     &#931; p&#8336; f(x&#8336;) = 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)---- make a constraint from function and constant-constraint :: Floating a => ExpectationFunction a -> a -> Constraint a-constraint = (,)---- The average constraint-average :: Floating a => a -> Constraint a-average m = constraint id m---- The variance constraint-variance :: Floating a => a -> Constraint a-variance sigma = constraint (^(2 :: Int)) sigma---- | Most general solver---   This will solve the langrangian by added the constraint that the ---   probabilities must add up to zero.---   This is the slowest but most flexible method. --generalMaxent :: (forall a. Floating 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 :: Floating a => [a] -> a-   obj = generalObjectiveFunc objInput--   values :: Floating a => [a]-   values = fst params--   constraints :: Floating a => [(GeneralConstraint a, a)]-   constraints = snd params--   fsmoments :: Floating a => ([GeneralConstraint a], [a])-   fsmoments = unzip constraints -   -   --The new constraint is always needed for probability problems-   -   sumToOne :: Floating a => GeneralConstraint a-   sumToOne _ xs = sumMap id xs-   -   gcons' :: Floating a => [GeneralConstraint a]-   gcons' = fst fsmoments-   -   gcons :: Floating a => [GeneralConstraint a]-   gcons = sumToOne : gcons'-   -   moments :: Floating a => [a]-   moments = 1 : snd fsmoments-   -   -   objInput :: Floating 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)-       -       ---- | This is for the linear case Ax = b ---   @x@ in this situation is the vector of probablities.---  ---  For example.--- --- @---   maxentLinear ([1,1,1], ([[0.85, 0.1, 0.05], [0.25, 0.5, 0.25], [0.05, 0.1, 0.85]], [0.29, 0.42, 0.29]))--- @------ Right [0.1, 0.8, 0.1]--- --- To be honest I am not sure why I can't use the 'maxent' version to solve--- this type of problem, but it doesn't work. I'm still learning--- -maxentLinear :: (forall a. Floating a => ([a], ([[a]], [a]))) -- ^ The values and a matrix A and column vector b-      -> Either (Result, Statistics) [Double] -- ^ Either the a discription of what wrong or the probability distribution -maxentLinear params = result where-   obj :: Floating a => [a] -> a-   obj = uncurry (objectiveFuncLinear values) fsmoments--   values :: Floating a => [a]-   values = fst params--   constraints :: Floating a => ([[a]], [a])-   constraints = snd params--   fsmoments :: Floating a => ([[a]], [a])-   fsmoments = constraints --   fs :: [[Double]]-   fs = fst fsmoments--   -- hmm maybe there is a better way to get rid of the defaulting-   guess = U.fromList $ replicate -       (length fs) (2.0 :: Double) --   result = case unsafePerformIO (optimize defaultParameters 0.00001 guess -                       (toFunction obj)-                       (toGradient obj)-                       Nothing) of-       (vs, ToleranceStatisfied, _) -> Right $ linProbs fs (S.toList vs)-       (_, x, y) -> Left (x, y)----- | The main entry point for computing discrete maximum entropy distributions.---   Where the constraints are all moment constraints. -maxent :: (forall a. Floating 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 :: 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 :: [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-        (vs, ToleranceStatisfied, _) -> Right $ probs values fs (S.toList vs)-        (_, x, y) -> Left (x, y)----    -    ---
+ src/Numeric/MaxEnt.hs view
@@ -0,0 +1,60 @@+-- |+-- The maximum entropy method, or MAXENT, is variational approach for computing probability +-- distributions given a list of moment, or expected value, constraints.+-- +-- Here are some links for background info.+-- A good overview of applications:+-- <http://cmm.cit.nih.gov/maxent/letsgo.html>+-- On the idea of maximum entropy in general: +-- <http://en.wikipedia.org/wiki/Principle_of_maximum_entropy>+--  +-- +-- Use this package to compute discrete maximum entropy distributions over a list of values and+-- list of constraints.+-- +-- 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]+-- +-- 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]+-- +-- 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.+-- +module Numeric.MaxEnt (+    Constraint,+    ExpectationFunction,+    constraint,+    average,+    variance,+    -- ** Classic moment based+    maxent,+    -- ** General+    GeneralConstraint,+    generalMaxent, +    -- ** Linear+    LinearConstraints(..),+    linear,+) where+import Numeric.MaxEnt.Internal (Constraint,+                        ExpectationFunction,+                        constraint,+                        average,+                        variance,+                        maxent,+                        generalMaxent,+                        GeneralConstraint)+import Numeric.MaxEnt.Linear (linear, LinearConstraints(..))+++++
+ src/Numeric/MaxEnt/ConjugateGradient.hs view
@@ -0,0 +1,50 @@+{-# LANGUAGE TupleSections, Rank2Types #-}+module Numeric.MaxEnt.ConjugateGradient 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)+++dot :: Num a => [a] -> [a] -> a+dot x y = sum . zipWith (*) x $ y++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+++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+      -> Int+      -> (forall a. RealFloat a => [a] -> a) +      -> Either (Result, Statistics) [Double]+solve percision 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+       (_, x, y) -> Left (x, y)
+ src/Numeric/MaxEnt/Internal.hs view
@@ -0,0 +1,210 @@+{-# 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+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 +--+-- @+--     &#931; p&#8336; f(x&#8336;) = 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)++++    +    +++
+ src/Numeric/MaxEnt/Linear.hs view
@@ -0,0 +1,115 @@+{-# LANGUAGE TupleSections, Rank2Types #-}+module Numeric.MaxEnt.Linear where+import Numeric.MaxEnt.ConjugateGradient+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++probs :: (RealFloat b) => [[b]] -> [b] -> [b]    +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 {+        matrix :: [[a]], +        output :: [a]+    }+    deriving (Show, Eq)++-- | This is for the linear case Ax = b +--   @x@ in this situation is the vector of probablities.+--  +--  Consider the 1 dimensional circular convolution using hmatrix.+--  +--  >>> import Numeric.LinearAlgebra+--  >>> fromLists [[0.68, 0.22, 0.1], [0.1, 0.68, 0.22], [0.22, 0.1, 0.68]] <> fromLists [[0.2], [0.5], [0.3]]+--  (3><1) [0.276, 0.426, 0.298]   +-- +--   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]+-- +--+-- +-- +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 tolerance constraints = result where+   obj :: RealFloat a => [a] -> a+   obj = objectiveFunc (matrix constraints) (output constraints) ++   as :: [[Double]]+   as = matrix constraints+   +   result = probs as <$> solve tolerance (length as) 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)))++  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+++++++-}++++++++++  +  ++++++
+ tests/Main.hs view
@@ -0,0 +1,18 @@+module Main where+import Test.Framework (defaultMain, testGroup, defaultMainWithArgs)+import Test.Framework.Providers.HUnit+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Numeric.MaxEnt+import LinearTesting++-- One test is that a list convolved with a +-- guassian and then infered with a gaussian should have less entropy then it+-- started with++main = defaultMain [+        testGroup "Linear Tests" [+            testProperty "solvableSystemsAreSolvable" solvableSystemsAreSolvable,+            testProperty "probsSumToOne" probsSumToOne,+            testProperty "solutionFitsConstraints" solutionFitsConstraints   +        ]+    ]