{-# LANGUAGE DeriveDataTypeable #-} {-| Module : Data.Number.ER.RnToRm.BisectionTree Description : hierarchical domain partitions Copyright : (c) 2007-2008 Michal Konecny License : BSD3 Maintainer : mik@konecny.aow.cz Stability : experimental Portability : portable Defines a representation for recursive bisections of @R^n@ by hyperplanes, each of which is perpendicular to a base axis. Arbitrary data can be associated with the sections of a partition. To be imported qualified, usually with prefix BISTR. -} module Data.Number.ER.RnToRm.BisectionTree ( BisectionTree(..), Depth, ValueSplitter, ValueCombiner, isLeaf, const, removeVars, sync2, syncMany, setDepth, split, mapWithDom, mapLeaves, doBistr, doMap, doMapLeaves, combineWith, collectValues, collectDomValues, compare, lookupSubtreeDoms ) where import Prelude hiding (const, map, compare) import qualified Prelude import qualified Data.Number.ER.Real.Approx as RA import qualified Data.Number.ER.BasicTypes.DomainBox as DBox import Data.Number.ER.BasicTypes.DomainBox (VariableID(..), DomainBox, DomainBoxMappable, DomainIntBox) import Data.Number.ER.BasicTypes import Data.Number.ER.Misc import Data.Number.ER.ShowHTML import qualified Text.Html as H import Data.Typeable import Data.Generics.Basics import Data.Binary --import BinaryDerive import Data.Maybe {-| * The root of the tree often represents the whole @R^n@. * Each node splits the parent's space into two using a specified variable (ie direction) and an optional splitting point. * By default, a split is taken at the point defined by the method 'RA.bisect'. -} data BisectionTree box varid d v = Leaf { bistrDepth :: Depth, bistrDom :: box, -- ^ domain bistrVal :: v -- ^ value estimate } | Node { bistrDepth :: Depth, -- ^ depth of this node bistrDom :: box, -- ^ domain bistrDir :: varid, -- ^ direction to split in bistrPt :: d, -- ^ point that the split is at bistrLO :: BisectionTree box varid d v, -- ^ the half towards -Infty in split dir bistrHI :: BisectionTree box varid d v -- ^ the half towards +Infty in split dir } deriving (Typeable, Data) type Depth = Int {- the following has been generated by BinaryDerive -} instance (Binary a, Binary b, Binary c, Binary d) => Binary (BisectionTree a b c d) where put (Leaf a b c) = putWord8 0 >> put a >> put b >> put c put (Node a b c d e f) = putWord8 1 >> put a >> put b >> put c >> put d >> put e >> put f get = do tag_ <- getWord8 case tag_ of 0 -> get >>= \a -> get >>= \b -> get >>= \c -> return (Leaf a b c) 1 -> get >>= \a -> get >>= \b -> get >>= \c -> get >>= \d -> get >>= \e -> get >>= \f -> return (Node a b c d e f) _ -> fail "no parse" {- the above has been generated by BinaryDerive -} instance (VariableID varid, Show d, Show v, DomainBox box varid d) => Show (BisectionTree box varid d v) where show = showBisectionTree show showBisectionTree showValue = showB where showB (Leaf depth domB val) = "\n" ++ (concat (replicate (depth * 2) ".")) ++ "o " ++ (concatWith "," (Prelude.map showVD $ DBox.toList domB)) ++ " |---> " ++ showValue val showB (Node depth domB dir pt lo hi) = "\n" ++ (concat (replicate (depth * 2) ".")) ++ "o " ++ (concatWith "," (Prelude.map showVD $ DBox.toList domB)) ++ " //" ++ showVar dir ++ "\\\\" ++ (concat $ Prelude.map (showBisectionTree showValue) [lo,hi]) showVD (v,d) = showVar v ++ "->" ++ show d instance (Show d, H.HTML v, DomainBox box varid d) => H.HTML (BisectionTree box varid d v) where toHtml (Leaf depth domB val) = H.toHtmlFromList $ [ H.toHtml $ concatWith "," (Prelude.map showVD $ DBox.toList domB) , H.primHtml " → " , H.toHtml val ] where showVD (v,d) = showVar v ++ " in " ++ show d toHtml (Node depth domB dir pt lo hi) = H.toHtml $ besidesTable [H.border 2] [ abovesTable [] [H.toHtml $ "(" ++ showVar dir ++ ")"] , abovesTable [] [H.toHtml lo, H.toHtml hi] ] isLeaf :: BisectionTree box varid d v -> Bool isLeaf (Leaf _ _ _) = True isLeaf (Node _ _ _ _ _ _) = False const :: -- (DomainIntBox box varid d) => box -> v -> BisectionTree box varid d v const domB value = Leaf 0 domB value {-| value splitter function - parameters are: depth, domain of value, value, variable to split by, point to split at; returns the two split values -} type ValueSplitter box varid d v = (EffortIndex -> Depth -> box -> v -> varid -> d -> (v,v)) type ValueCombiner box varid d v = (EffortIndex -> Depth -> (BisectionTree box varid d v) -> v) setDepth :: Depth -> BisectionTree box varid d v -> BisectionTree box varid d v setDepth depth bistr | isLeaf bistr = bistr { bistrDepth = depth } | otherwise = bistr { bistrLO = setDepth depthInc $ bistrLO bistr, bistrHI = setDepth depthInc $ bistrHI bistr } where depthInc = depth + 1 split :: (RA.ERIntApprox d, DomainBox box varid d) => ValueSplitter box varid d v -> EffortIndex -> varid {-^ variable @x@ to split by -} -> d {-^ point in domain of @x@ to split at -} -> box {-^ domain to lookup @x@ in if tree's domain does not have @x@ -} -> BisectionTree box varid d v -> BisectionTree box varid d v split valSplitter ix splitDir splitPt fallbackDom bistr = resultBistr where resultBistr = spl bistr spl (Leaf depth domB val) = Node depth domB splitDir splitPt childLO childHI where childLO = Leaf depthInc domLO valLO childHI = Leaf depthInc domHI valHI (valLO, valHI) = valSplitter ix depth domB val splitDir splitPt depthInc = depth + 1 domLO = DBox.insert splitDir dirDomLO domB domHI = DBox.insert splitDir dirDomHI domB (dirDomLO, dirDomHI) = RA.bisectDomain (Just splitPt) dirDom dirDom = DBox.findWithDefault (DBox.lookup "BisectionTree: split: fallbackDom: " splitDir fallbackDom) splitDir domB spl bistr@(Node depth domB dir pt childLO childHI) | dir == splitDir = case RA.compareReals pt splitPt of Just LT -> -- split on lower half Node depth domB dir pt (Node depthInc domChildLO splitDir splitPt childLOsplitLO childLOsplitHI) childHI Just GT -> -- split on higher half Node depth domB dir pt childLO (Node depthInc domChildHI splitDir splitPt childHIsplitLO childHIsplitHI) _ -> bistr | otherwise = -- splitDir < dir = Node depth domB dir pt (Node depthInc domChildLO splitDir splitPt childLOsplitLO childLOsplitHI) (Node depthInc domChildHI splitDir splitPt childHIsplitLO childHIsplitHI) -- | dir < splitDir = -- Node depth domB dir childLOsplit childHIsplit where depthInc = depth + 1 domChildLO = bistrDom childLO domChildHI = bistrDom childHI childLOsplit@(Node _ _ _ _ childLOsplitLO childLOsplitHI) = spl childLO childHIsplit@(Node _ _ _ _ childHIsplitLO childHIsplitHI) = spl childHI {-| Apply a function to all values, thus creating a new tree. -} mapWithDom :: (DomainBox box varid d) => (box -> v1 -> v2) -> BisectionTree box varid d v1 -> BisectionTree box varid d v2 mapWithDom f bistr@(Leaf _ domB val) = bistr { bistrVal = f domB val } mapWithDom f bistr@(Node _ _ _ _ cLO cHI) = bistr { bistrLO = mapWithDom f cLO, bistrHI = mapWithDom f cHI } {-| Apply a function to all values, thus creating a new tree. -} mapLeaves :: (BisectionTree box varid d v1 -> BisectionTree box varid d v2) -> BisectionTree box varid d v1 -> BisectionTree box varid d v2 mapLeaves f bistr@(Leaf _ domB val) = f bistr mapLeaves f bistr@(Node _ _ _ _ cLO cHI) = bistr { bistrLO = mapLeaves f cLO, bistrHI = mapLeaves f cHI } {-| Apply a function to all values, thus creating a list of new trees. -} mapMultiLeaves :: (BisectionTree box varid d v1 -> [BisectionTree box varid d v2]) -> BisectionTree box varid d v1 -> [BisectionTree box varid d v2] mapMultiLeaves f bistr@(Leaf _ domB val) = f bistr mapMultiLeaves f bistr@(Node _ _ _ _ cLO cHI) = Prelude.map (replaceChildren bistr) $ zip (mapMultiLeaves f cLO) (mapMultiLeaves f cHI) where replaceChildren bistr (newLO, newHI) = bistr { bistrLO = newLO, bistrHI = newHI } {-| Perform a given action on all branches of a bisection tree, left to right. (optionally now going below the given depth) -} doBistr :: (box -> v -> IO ()) -> Maybe Int -> BisectionTree box varid d v -> IO () doBistr f Nothing bistr = m bistr where m (Node _ _ _ _ lo hi) = do m lo m hi m (Leaf _ domB val) = f domB val doBistr f (Just maxDepth) bistr = m maxDepth bistr where m maxDepth (Node depth domB _ _ lo hi) | maxDepth > 0 = do m (maxDepth - 1) lo m (maxDepth - 1) hi | otherwise = error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached" -- m err (Leaf depth domB val) -- where -- val = head $ collectValues lo -- err = m _ (Leaf _ domB val) = f domB val {-| Perform a given action on all branches of a bisection tree, left to right. (optionally now going below the given depth) -} doMap :: (Depth -> box -> v -> IO v) -> Maybe Int -> BisectionTree box varid d v -> IO (BisectionTree box varid d v) doMap f Nothing bistr = m bistr where m bistr@(Node _ _ _ _ lo hi) = do newLo <- m lo newHi <- m hi return $ bistr { bistrLO = newLo, bistrHI = newHi } m bistr@(Leaf depth domB val) = do newVal <- f depth domB val return $ bistr { bistrVal = newVal } doMap f (Just maxDepth) bistr = m maxDepth bistr where m maxDepth bistr@(Node depth domB _ _ lo hi) | maxDepth > 0 = do newLo <- m (maxDepth - 1) lo newHi <- m (maxDepth - 1) hi return $ bistr { bistrLO = newLo, bistrHI = newHi } | otherwise = error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached" -- m err (Leaf depth domB val) -- where -- val = head $ collectValues lo -- err = m _ bistr@(Leaf depth domB val) = do newVal <- f depth domB val return $ bistr { bistrVal = newVal } {-| Perform a given action on all branches of a bisection tree, left to right with the option of further branching the tree. (optionally now going below the given depth) -} doMapLeaves :: (BisectionTree box varid d v -> IO (BisectionTree box varid d v)) -> Maybe Int -> BisectionTree box varid d v -> IO (BisectionTree box varid d v) doMapLeaves f Nothing bistr = m bistr where m bistr@(Node _ _ _ _ lo hi) = do newLo <- m lo newHi <- m hi return $ bistr { bistrLO = newLo, bistrHI = newHi } m bistr@(Leaf depth domB val) = do f bistr doMapLeaves f (Just maxDepth) bistr = m maxDepth bistr where m maxDepth bistr@(Node depth domB _ _ lo hi) | maxDepth > 0 = do newLo <- m (maxDepth - 1) lo newHi <- m (maxDepth - 1) hi return $ bistr { bistrLO = newLo, bistrHI = newHi } | otherwise = error $ "BisectionTree: doBistr: maxDepth (=" ++ show maxDepth ++ ") breached" -- m err (Leaf depth domB val) -- where -- val = head $ collectValues lo -- err = m _ bistr@(Leaf depth domB val) = do f bistr removeVars :: (RA.ERIntApprox d, DomainIntBox box varid d, DomainBoxMappable box box varid d d) => box -> BisectionTree box varid d v -> BisectionTree box varid d v removeVars substitutions bistr = aux (bistrDepth bistr) bistr where aux depth (Leaf _ domB val) = Leaf depth domNoVars val where domNoVars = DBox.difference domB substitutions aux depth (Node _ domB v pt lo hi) | v `DBox.member` substitutions = case (vVal `RA.refines` vDomLO, vVal `RA.refines` vDomHI) of (True, _) -> aux depth lo (_, True) -> aux depth hi | otherwise = Node depth domNoVars v pt loNoVars hiNoVars where vVal = DBox.lookup loc v substitutions vDomLO = DBox.lookup loc v $ bistrDom lo vDomHI = DBox.lookup loc v $ bistrDom hi loc = "RnToRm.BisectionTree: removeVars: " domNoVars = DBox.difference domB substitutions loNoVars = aux (depth + 1) lo hiNoVars = aux (depth + 1) hi {-| Ensure both trees have equal structure at the top level: either they are all leaves or they all split at the same direction with the same splitting point. Also, unify the domains at the top level. -} sync2 :: (RA.ERIntApprox d, DomainIntBox box varid d) => ValueSplitter box varid d v1 -> ValueSplitter box varid d v2 -> EffortIndex -> BisectionTree box varid d v1 -> BisectionTree box varid d v2 -> (BisectionTree box varid d v1, BisectionTree box varid d v2) sync2 valSplitter1 valSplitter2 ix bistr1 bistr2 = case getPt bistr1 bistr2 of Nothing -> unifyDom bistr1 bistr2 Just (var, pt, domB) -> unifyDom (split valSplitter1 ix var pt domB bistr1) (split valSplitter2 ix var pt domB bistr2) where getPt bistr1 bistr2 | isLeaf bistr1 && isLeaf bistr2 = Nothing | isLeaf bistr1 = Just (bistrDir bistr2, bistrPt bistr2, bistrDom bistr2) | otherwise = Just (bistrDir bistr1, bistrPt bistr1, bistrDom bistr1) unifyDom bistr1 bistr2 = (bistr1 { bistrDom = domB }, bistr2 { bistrDom = domB }) where domB = DBox.unify "RnToRm.BisectionTree: sync: " dom1 dom2 dom1 = bistrDom bistr1 dom2 = bistrDom bistr2 {-| Ensure all the trees have equal structure at the top level: either they are all leaves or they all split at the same direction with the same splitting point. Also, unify the domains at the top level. -} syncMany :: (RA.ERIntApprox d, DomainIntBox box varid d) => ValueSplitter box varid d v -> EffortIndex -> [BisectionTree box varid d v] -> [BisectionTree box varid d v] syncMany valSplitter ix bistrs = case getPt bistrs of Nothing -> unifyDom bistrs Just (var, pt, domB) -> unifyDom $ Prelude.map (split valSplitter ix var pt domB) bistrs where getPt [] = Nothing getPt (bistr : rest) | isLeaf bistr = getPt rest | otherwise = Just (bistrDir bistr, bistrPt bistr, bistrDom bistr) unifyDom bistrs = Prelude.map (setDom domB) bistrs where setDom domB bistr = bistr { bistrDom = domB } domB = foldl (DBox.unify "RnToRm.BisectionTree: sync: ") DBox.noinfo $ Prelude.map bistrDom bistrs {-| Combine two bisection trees using a given value combining function. Where necessary, leaves are split so that the resulting tree's structure is the union of the two argument tree structures. Such splitting of values in leaves is performed by the provided functions. -} combineWith :: (RA.ERIntApprox d, DomainIntBox box varid d) => ValueSplitter box varid d v1 {-^ value splitter function for tree 1 -} -> ValueSplitter box varid d v2 {-^ value splitter function for tree 2 -} -> (box -> v1 -> v2 -> (Maybe res, aux)) {-^ partial function to combine values with -} -> EffortIndex -> (BisectionTree box varid d v1) -> (BisectionTree box varid d v2) -> (Maybe (BisectionTree box varid d res), [aux]) combineWith valSplitter1 valSplitter2 f ix bistr1 bistr2 = combineAux bistr1sync bistr2sync where (bistr1sync, bistr2sync) = sync2 valSplitter1 valSplitter2 ix bistr1 bistr2 combineAux bistr1@(Leaf _ domB val1) bistr2@(Leaf _ _ val2) = case f domB val1 val2 of (Nothing, aux) -> (Nothing, [aux]) (Just val, aux) -> (Just $ bistr1 { bistrVal = val }, [aux]) combineAux bistr1@(Node _ domB _ _ lo1 hi1) bistr2@(Node _ _ _ _ lo2 hi2) = ( Just $ bistr1 { bistrLO = fromJust mbistrLO, bistrHI = fromJust mbistrHI } , auxLO ++ auxHI ) where (mbistrLO, auxLO) = combineAux lo1Sync lo2Sync (mbistrHI, auxHI) = combineAux hi1Sync hi2Sync (lo1Sync, lo2Sync) = sync2 valSplitter1 valSplitter2 ix lo1 lo2 (hi1Sync, hi2Sync) = sync2 valSplitter1 valSplitter2 ix hi1 hi2 {-| return all values in leafs (except those within some CE subtree) as a list (from the leftmost to the rightmost) -} collectValues :: BisectionTree box varid b a -> [a] collectValues (Leaf _ _ val) = [val] collectValues (Node _ _ _ _ cLO cHI) = (collectValues cLO) ++ (collectValues cHI) {-| return all values in leafs (except those within some CE subtree) as a list (from the leftmost to the rightmost) -} collectDomValues :: BisectionTree box varid d v -> [(box, v)] collectDomValues (Leaf _ domB val) = [(domB,val)] collectDomValues (Node _ _ _ _ cLO cHI) = (collectDomValues cLO) ++ (collectDomValues cHI) {-| linear ordering on bisection trees -} compare :: (Ord varid, DomainBox box varid d) => (d -> d -> Ordering) -> (v -> v -> Ordering) -> (BisectionTree box varid d v) -> (BisectionTree box varid d v) -> Ordering compare compDoms compValues (Leaf _ _ _) (Node _ _ _ _ _ _) = LT compare compDoms compValues (Node _ _ _ _ _ _) (Leaf _ _ _) = GT compare compDoms compValues (Leaf depth1 dom1 val1) (Leaf depth2 dom2 val2) = compareComposeMany [ Prelude.compare depth1 depth2, DBox.compare compDoms dom1 dom2, compValues val1 val2 ] compare compDoms compValues (Node depth1 dom1 dir1 pt1 lo1 hi1) (Node depth2 dom2 dir2 pt2 lo2 hi2) = compareComposeMany [ Prelude.compare dir1 dir2, compDoms pt1 pt2, compare compDoms compValues lo1 lo2, compare compDoms compValues hi1 hi2 ] {-| lookup all maximal subtrees whose domain intersect the given rectangle -} lookupSubtreeDoms :: (RA.ERIntApprox d, DomainBox box varid d) => (BisectionTree box varid d v) -> box {-^ domain to look up within the tree -} -> [BisectionTree box varid d v] lookupSubtreeDoms origBistr domB = lk origBistr where lk bistr@(Leaf _ _ _) = [bistr] lk bistr@(Node _ _ _ _ lo hi) | loDisjoint = lk hi | hiDisjoint = lk lo | otherwise = (lk lo) ++ (lk hi) where loDisjoint = and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.isDisjoint) domB domLO hiDisjoint = and $ Prelude.map snd $ DBox.zipWithDefault RA.bottomApprox (RA.isDisjoint) domB domHI domLO = bistrDom lo domHI = bistrDom hi {-| Update a value on a given sub-domain, bisecting the tree if necessary to obtain a better fit for the domain, but not below a given depth limit. With multiple domain dimensions, split the domain according to `DBox.bestSplit'. -} updateVal :: (RA.ERIntApprox d, DomainIntBox box varid d, DomainBoxMappable box box varid d d) => ValueSplitter box varid d v -> EffortIndex -> Depth {-^ depth limit -} -> box {-^ domain to update on -} -> (box -> v -> v) {-^ how to update values that intersect the above domain -} -> (BisectionTree box varid d v) -> (BisectionTree box varid d v) updateVal valSplitter ix maxDepth updateDom updateFn bistr = upd bistr where upd bistr | noOverlap = bistr | edgeTouch && (isLeaf bistr) = updateLeaf bistr -- assuming we can update values on edges without -- influence on the interior | insideUpdateDom = mapLeaves updateLeaf bistr | depth >= maxDepth = mapLeaves updateLeaf bistr | otherwise = -- divide and conquer: Node depth domB dir pt bistrLdone bistrRdone where updateLeaf bistr = bistr { bistrVal = updateFn (bistrDom bistr) (bistrVal bistr) } noOverlap = or $ Prelude.map (not . RA.isConsistent) $ DBox.elems domOverlap domOverlap = DBox.intersectionWith (RA./\) domB updateDom insideUpdateDom = and $ Prelude.map snd $ DBox.zipWith RA.refines domB updateDom edgeTouch = and $ Prelude.map snd $ DBox.zipWithDefaultSecond RA.bottomApprox endPointTouch domB updateDom endPointTouch i1 i2 = i1L == i2R || i1R == i2L where (==) = RA.eqSingletons (i1L, i1R) = RA.bounds i1 (i2L, i2R) = RA.bounds i2 depth = bistrDepth bistr domB = bistrDom bistr bistrLdone = upd bistrL bistrRdone = upd bistrR (Node _ _ _ _ bistrL bistrR) | (isLeaf bistr) = split valSplitter ix dir pt DBox.noinfo bistr | otherwise = bistr (dir, (_,pt)) = DBox.bestSplit domB