module Bayes.Factor(
Factor(..)
, isomorphicFactor
, normedFactor
, Set(..)
, BayesianDiscreteVariable(..)
, Vertex(..)
, DV
, DVI
, setDVValue
, instantiationValue
, instantiationVariable
, variableVertex
, (=:)
, forAllInstantiations
, factorFromInstantiation
, changeVariableOrder
, CPT
, testProductProject_prop
, testScale_prop
, testProjectCommut_prop
, testScalarProduct_prop
, testProjectionToScalar_prop
) where
import qualified Data.Vector.Unboxed as V
import Data.Vector.Unboxed((!))
import Data.Maybe(fromJust,mapMaybe,isJust)
import qualified Data.List as L
import Text.PrettyPrint.Boxes hiding((//))
import Test.QuickCheck
import Test.QuickCheck.Arbitrary
import qualified Data.IntMap as IM
import Control.Monad
import System.Random(Random)
import Data.List(partition)
import Bayes.PrivateTypes
class LabeledVertex l where
variableVertex :: l -> Vertex
factorFromInstantiation :: Factor f => DVI Int -> f
factorFromInstantiation (DVI dv a) =
let setValue i = if i == a then 1.0 else 0.0
in
fromJust . factorWithVariables [dv] . map (setValue) $ [0..dimension dv1]
instance LabeledVertex (DVI a) where
variableVertex (DVI v _) = variableVertex v
instance LabeledVertex DV where
variableVertex (DV v _) = v
normedFactor :: Factor f => f -> f
normedFactor f = factorDivide f (factorNorm f)
class Factor f where
isScalarFactor :: f -> Bool
emptyFactor :: f
containsVariable :: f -> DV -> Bool
factorVariables :: f -> [DV]
factorMainVariable :: f -> DV
factorMainVariable f = let vars = factorVariables f
in
case vars of
[] -> error "Can't get the main variable of a scalar factor"
(h:_) -> h
factorWithVariables :: [DV] -> [Double] -> Maybe f
factorValue :: f -> [DVI Int] -> Double
variablePosition :: f -> DV -> Maybe Int
factorDimension :: f -> Int
factorNorm :: f -> Double
factorScale :: Double -> f -> f
factorFromScalar :: Double -> f
evidenceFrom :: [DVI Int] -> Maybe f
factorDivide :: f -> Double -> f
factorDivide f d = (1.0 / d) `factorScale` f
factorToList :: f -> [Double]
factorProduct :: [f] -> f
factorProjectOut :: [DV] -> f -> f
factorProjectTo :: [DV] -> f -> f
factorProjectTo s f =
let alls = factorVariables f
s' = alls `difference` s
in
factorProjectOut s' f
changeVariableOrder :: DVSet s
-> DVSet s'
-> [Double]
-> [Double]
changeVariableOrder (DVSet oldOrder) newOrder oldValues =
let oldFactor = fromJust $ factorWithVariables oldOrder oldValues :: CPT
in
[factorValue oldFactor i | i <- forAllInstantiations newOrder]
data CPT = CPT {
dimensions :: ![DV]
, mapping :: !(IM.IntMap Int)
, values :: !(V.Vector Double)
}
| Scalar !Double
debugCPT (Scalar d) = do
putStrLn "SCALAR CPT"
print d
putStrLn ""
debugCPT (CPT d m v) = do
putStrLn "CPT"
print d
putStrLn ""
print m
putStrLn ""
print v
putStrLn ""
quickCheckVertexSize :: Int -> Int
quickCheckVertexSize 0 = 2
quickCheckVertexSize 1 = 2
quickCheckVertexSize 2 = 2
quickCheckVertexSize _ = 2
whileIn :: (Arbitrary a, Eq a) => [a] -> Gen a -> Gen a
whileIn l m = do
newVal <- m
if newVal `elem` l
then
whileIn l m
else
return newVal
generateWithoutReplacement :: (Random a, Arbitrary a, Eq a)
=> Int
-> (a,a)
-> Gen [a]
generateWithoutReplacement n b | n == 1 = generateSingle b
| n > 1 = generateMultiple n b
| otherwise = return []
where
generateSingle b = do
r <- choose b
return [r]
generateMultiple n b = do
l <- generateWithoutReplacement (n1) b
newelem <- whileIn l $ choose b
return (newelem:l)
instance Arbitrary CPT where
arbitrary = do
nbVertex <- choose (1,4) :: Gen Int
vertexNumbers <- generateWithoutReplacement nbVertex (0,50)
let dimensions = map (\i -> (DV (Vertex i) (quickCheckVertexSize i))) vertexNumbers
let valuelen = product (map dimension dimensions)
rndValues <- vectorOf valuelen (choose (0.0,1.0) :: Gen Double)
return . fromJust . factorWithVariables dimensions $ rndValues
nearlyEqual :: Double -> Double -> Bool
nearlyEqual a b =
let absA = abs a
absB = abs b
diff = abs (ab)
epsilon = 2e-5
in
case (a,b) of
(x,y) | x == y -> True
| x*y == 0 -> diff < (epsilon * epsilon)
| otherwise -> diff / (absA + absB) < epsilon
testScale_prop :: Double -> CPT -> Bool
testScale_prop s f = (factorNorm (s `factorScale` f)) `nearlyEqual` (s * (factorNorm f))
testProductProject_prop :: CPT -> CPT -> Property
testProductProject_prop fa fb = isEmpty ((factorVariables fa) `intersection` (factorVariables fb)) ==>
let r = factorProjectOut (factorVariables fb) (factorProduct [fa,fb])
fa' = r `factorDivide` (factorNorm fb)
in
fa' `isomorphicFactor` fa
testScalarProduct_prop :: Double -> CPT -> Bool
testScalarProduct_prop v f = (factorProduct [(Scalar v),f]) `isomorphicFactor` (v `factorScale` f)
testProjectionToScalar_prop :: CPT -> Bool
testProjectionToScalar_prop f =
let allVars = factorVariables f
in
(factorProjectOut allVars f) `isomorphicFactor` (factorFromScalar (factorNorm f))
testProjectCommut_prop:: CPT -> Property
testProjectCommut_prop f = length (factorVariables f) >= 3 ==>
let a = take 1 (factorVariables f)
b = take 1 . drop 1 $ factorVariables f
commuta = factorProjectOut a (factorProjectOut b f)
commutb = factorProjectOut b (factorProjectOut a f)
in
commuta `isomorphicFactor` commutb
isomorphicFactor :: Factor f => f -> f -> Bool
isomorphicFactor fa fb = maybe False (const True) $ do
let sa = factorVariables fa
sb = factorVariables fb
va = DVSet sa
vb = DVSet sb
guard (sa `equal` sb)
guard (factorDimension fa == factorDimension fb)
guard $ and [factorValue fa ia `nearlyEqual` factorValue fb ia | ia <- forAllInstantiations va]
return ()
vname :: Int -> DVI Int -> Box
vname vc i = text $ "v" ++ show vc ++ "=" ++ show (instantiationValue i)
dispFactor :: Factor f => f -> DV -> [DVI Int] -> [DV] -> Box
dispFactor cpt h c [] =
let dstIndexes = allInstantiationsForOneVariable h
dependentIndexes = reverse c
factorValueAtPosition p =
let v = factorValue cpt p
in
text . show $ v
in
vsep 0 center1 . map (factorValueAtPosition . (:dependentIndexes)) $ dstIndexes
dispFactor cpt dst c (h@(DV (Vertex vc) i):l) =
let allInst = allInstantiationsForOneVariable h
in
hsep 1 top . map (\i -> vcat center1 [vname vc i,dispFactor cpt dst (i:c) l]) $ allInst
instance Show CPT where
show (Scalar v) = "\nScalar Factor:\n" ++ show v
show c@(CPT [] _ v) = "\nEmpty CPT:\n"
show c@(CPT d _ v) =
let h@(DV (Vertex vc) _) = head d
table = dispFactor c h [] (tail d)
dstIndexes = map head (forAllInstantiations . DVSet $ [h])
dstColumn = vcat center1 $ replicate (length d 1) (text "") ++ map (vname vc) dstIndexes
in
"\n" ++ show d ++ "\n" ++ render (hsep 1 top [dstColumn,table])
instance Factor CPT where
factorToList (Scalar v) = [v]
factorToList (CPT _ _ v) = V.toList v
emptyFactor = emptyCPT
isScalarFactor (Scalar _) = True
isScalarFactor _ = False
factorFromScalar v = Scalar v
factorDimension f@(CPT _ _ _) = product . map dimension . factorVariables$ f
factorDimension _ = 1
containsVariable (CPT _ m _) (DV (Vertex i) _) = IM.member i m
containsVariable (Scalar _) _ = False
factorWithVariables = createCPTWithDims
factorVariables (CPT v _ _) = v
factorVariables (Scalar _) = []
factorNorm f@(CPT d _ vals) =
let vars = DVSet d
strides = indexStrides vars
in
sum [ vals!(indexPosition strides x) | x <- indicesForDomain vars]
factorNorm (Scalar v) = v
variablePosition (CPT _ m _) (DV (Vertex i) _) = IM.lookup i m
variablePosition (Scalar _) _ = Nothing
factorScale s (Scalar v) = Scalar (s*v)
factorScale s f@(CPT d _ vals) =
let vars = DVSet d
strides = indexStrides vars
newValues = map (s *) [ vals!(indexPosition strides x) | x <- indicesForDomain vars]
in
fromJust $ factorWithVariables (factorVariables f) newValues
factorValue (Scalar v) _ = v
factorValue f@(CPT d _ v) i =
let vars = DVSet d
(dv,pos) = instantiationDetails i
strides = indexStridesFor vars dv
in
v!(indexPosition strides pos)
evidenceFrom [] = Nothing
evidenceFrom l =
let (variables,index) = instantiationDetails l
DVSet nakedVars = variables
setValueForIndex i = if i == index then 1.0 else 0.0
in
factorWithVariables nakedVars . map setValueForIndex $ indicesForDomain variables
factorProduct [] = factorFromScalar 1.0
factorProduct l =
let allVars = DVSet $ L.foldl1' union . map factorVariables $ l
DVSet nakedVars = allVars
(scalars,cpts) = partition isScalarFactor l
stridesFromCPT (CPT d _ _) = indexStridesFor (DVSet d) allVars
ps = product . map (flip factorValue undefined) $ scalars
cptsStrides = map stridesFromCPT cpts
in
if L.null cpts
then
factorFromScalar ps
else
let getFactorValueAtIndex i (strides,factor@(CPT _ _ vals)) = vals!(indexPosition strides i)
instantiationProduct instantiation = product . map (getFactorValueAtIndex instantiation) $ (zip cptsStrides cpts)
values = [ps * instantiationProduct x | x <- indicesForDomain allVars]
in
values `seq` fromJust $ factorWithVariables nakedVars values
factorProjectOut _ f@(Scalar v) = f
factorProjectOut s f@(CPT d _ v) =
let variablesToSum = s
variablesToKeep = d `difference` s
keepSet = DVSet variablesToKeep
sumSet = DVSet variablesToSum
strides = indexStridesFor (DVSet d) (DVSet $ variablesToKeep ++ variablesToSum)
values = do
keepIndex <- indicesForDomain keepSet
let l = do
sumIndex <- indicesForDomain sumSet
return $ v!(indexPosition strides $ combineIndex strides keepIndex sumIndex)
return (sum l)
in
values `seq` fromJust $ factorWithVariables variablesToKeep values
combineIndex :: Strides s'' -> [Index s] -> [Index s'] -> [Index s'']
combineIndex _ la lb = map (Index . fromIndex) la ++ map (Index .fromIndex) lb
emptyCPT :: CPT
emptyCPT = CPT [] IM.empty V.empty
newtype Strides s = Strides [Int] deriving(Eq,Show)
indexStridesFor :: DVSet s
-> DVSet s'
-> Strides s'
indexStridesFor dr@(DVSet drvars) di@(DVSet divars) =
let Strides originStrides = indexStrides dr
reference = zip drvars originStrides
getNewStrides dv = maybe 0 id (lookup dv reference)
in
Strides $ map getNewStrides divars
indexStrides :: DVSet s -> Strides s
indexStrides d@(DVSet dvars) =
let dim = map dimension dvars
pos' = scanr (*) (1::Int) (tail dim)
in
Strides pos'
indexPosition :: Strides s -> [Index s] -> Int
indexPositions _ [] = 0
indexPosition (Strides pos') pos = sum . map (\(x,y) -> x * fromIndex y) $ (zip pos' pos)
createCPTWithDims :: [DV] -> [Double] -> Maybe CPT
createCPTWithDims dims values =
let createDVIndex i (DV (Vertex v) _) = (v,i)
m = IM.fromList . zipWith createDVIndex ([0,1..]::[Int]) $ dims
p = product (map dimension dims)
in
if length values == p
then
Just $ CPT dims m (V.fromList values)
else
Nothing