Cleans up main by moving some algorithms to Preorder.hs

app/Main.hs

@@ -6,9 +6,9 @@ import Mealy
 import MealyRefine
 import Merger
 import Partition
+import Preorder
 
-import Control.Monad (forM_, when)
-import Control.Monad.Trans.State.Strict
+import Control.Monad (forM_)
 import Data.Bifunctor
 import Data.List (sort, sortOn, intercalate)
 import Data.List.Ordered (nubSort)
@@ -76,31 +76,13 @@ main = do
     printPartition partition
     )
 
-  -- First we check eqiuvalent partitions, so that we only work on one
-  -- item in each equivalence class. This could be merged with the next
-  -- phase of checking refinement, and that would be faster. But this is
-  -- simpler.
-  let checkRelsFor o1 p1 =
-        forM_ projections (\(o2, p2) -> do
-          (repr, ls) <- get
-          -- We skip if o2 is equivelent to an earlier o
-          when (o1 < o2 && o2 `Map.notMember` repr) $ do
-            case isEquivalent p1 p2 of
-              -- Equivalent: let o2 point to o1
-              True -> put (Map.insert o2 o1 repr, ls)
-              False -> return ()
-          )
-      checkAllRels projs =
-        forM_ projs (\(o1, p1) -> do
-          -- First we check if o1 is equivalent to an earlier o
-          -- If so, we skip it. Else, we add it to the unique
-          -- components and compare to all others.
-          (repr, ls) <- get
-          when (o1 `Map.notMember` repr) $ do
-            put (repr, (o1, p1):ls)
-            checkRelsFor o1 p1
-          )
-      (equiv, uniqPartitions) = execState (checkAllRels projections) (Map.empty, [])
+  -- First we check for equivalent partitions, so that we skip redundant work.
+  let preord p1 p2 = toPreorder (comparePartitions p1 p2)
+      (equiv, uniqPartitions) = equivalenceClasses preord projections
+
+  putStrLn ""
+  putStrLn "Representatives"
+  print . fmap fst $ uniqPartitions
 
   putStrLn ""
   putStrLn "Equivalences"
@@ -108,35 +90,9 @@ main = do
     putStrLn $ "  " <> (show o2) <> " == " <> (show o1)
     )
 
-  -- Then we compare each pair of partitions. If one is a coarsening of
-  -- another, we can skip it later on. That is to say: we only want the
-  -- finest partitions.
-  let compareAll partitions =
-        forM_ partitions (\(o1, b1) ->
-          forM_ partitions (\(o2, b2) ->
-            when (o1 < o2) $ do
-              ls <- get
-              case comparePartitions b1 b2 of
-                Equivalent -> error "cannot happen"
-                Refinement -> put $ (o1, o2):ls
-                Coarsening -> put $ (o2, o1):ls
-                Incomparable -> return ()
-            )
-          )
-      rel = execState (compareAll uniqPartitions) []
-
-  putStrLn ""
-  putStrLn "Relation, coarser points to finer (bigger)"
-  forM_ rel (\(o1, o2) -> do
-    putStrLn $ "  " <> (show o2) <> " -> " <> (show o1)
-    )
-
-  -- Get rid of the coarser partitions
-  let lowElements = Set.fromList . fmap snd $ rel
-      allElements = Set.fromList . fmap fst $ uniqPartitions
-      topElements = Set.difference allElements lowElements
-      mods = Map.fromList uniqPartitions -- could be a lazy map
-      topMods = Map.assocs $ Map.restrictKeys mods topElements
+  -- Then we compare each pair of partitions. We only keep the finest
+  -- partitions, since the coarse ones don't provide value to us.
+  let (topMods, downSets) = maximalElements preord uniqPartitions
       foo (a, b) = (numBlocks b, a)
 
   putStrLn ""
@@ -145,21 +101,23 @@ main = do
     putStrLn $ "  " <> (show o) <> " has size " <> (show b)
     )
 
+  -- Then we try to combine paritions, so that we don't end up with
+  -- too many components. (Which would be too big to be useful.)
   let strategy MergerStats{..}
         | numberOfComponents <= 4 = Stop
         | otherwise = Continue
 
   projmap <- heuristicMerger topMods strategy
 
+  -- Now we are going to output the components we found.
   let equivInv = converseRelation equiv
-      relMap = Map.fromListWith (++) . fmap (second pure) $ rel
-      projmapN = zip projmap [1..]
+      projmapN = zip projmap [1 :: Int ..]
 
   forM_ projmapN (\((os, p), i) -> do
     let name = intercalate "x" os
         filename = "component" <> show i <> ".dot"
 
-        osWithRel = concat $ os:[Map.findWithDefault [] o relMap | o <- os]
+        osWithRel = concat $ os:[Map.findWithDefault [] o downSets | o <- os]
         osWithRelAndEquiv = concat $ osWithRel:[Map.findWithDefault [] o equivInv | o <- osWithRel]
         componentOutputs = Set.fromList osWithRelAndEquiv
         proj = projectToComponent (flip Set.member componentOutputs) machine
@@ -183,8 +141,18 @@ main = do
         initialBlock = state2block initialState
         -- Sorting on "/= initialBlock" puts the initialBlock in front
         initialFirst = sortOn (\(s,_,_,_) -> s /= initialBlock) transitionsBlocks
+
+        -- So far so good, `initialFirst` could serve as our output
+        -- But we do one more optimisation on the machine
+        -- We remove inputs, on which the machine does nothing
+        deadInputs0 = Map.fromListWith (++) . fmap (\(s,i,o,t) -> (i, [(s,o,t)])) $ initialFirst
+        deadInputs = Map.keysSet . Map.filter (all (\(s,o,t) -> s == t && o == Nothing)) $ deadInputs0
+        result = filter (\(_,i,_,_) -> i `Set.notMember` deadInputs) initialFirst
+
         -- Convert to a file
-        content = toString . mealyToDot name $ initialFirst
+        content = toString . mealyToDot name $ result
+
+    putStrLn $ "Dead inputs = " <> show (Set.size deadInputs)
 
     writeFile filename content
     )

mealy-decompose.cabal

@@ -30,7 +30,8 @@ library
     Mealy,
     MealyRefine,
     Merger,
-    Partition
+    Partition,
+    Preorder
   build-depends:
     vector
 

src/Partition.hs

@@ -3,14 +3,18 @@ module Partition
   , module Data.Partition
   ) where
 
+import Preorder
+
 import Control.Monad.Trans.State.Strict (runState, get, put)
+import Data.Coerce (coerce)
 import Data.Map.Strict qualified as Map
 import Data.Partition (Partition(..), numStates, blockOfState)
+import Data.Partition.Common (Block(..))
 import Data.Vector qualified as V
-import Unsafe.Coerce (unsafeCoerce)
 
--- Returns the coarsest partition which is finer than either input
--- i.e., the greatest lower bound
+-- | Returns the common refinement of two partitions. This is the coarsest
+-- partition which is finer than either input, i.e., the greatest lower bound.
+-- (If we put the finest partition on the top, and the coarsest on the bottom.)
 commonRefinement :: Partition -> Partition -> Partition
 commonRefinement p1 p2 =
   let n = numStates p1
@@ -25,7 +29,7 @@ commonRefinement p1 p2 =
             put (Map.insert key b m, succ b)
             return b
       (vect, (_, nextBlock)) = runState (V.generateM n blockAtIdx) (Map.empty, 0)
-  in Partition { numBlocks = unsafeCoerce nextBlock, stateAssignment = vect }
+  in Partition { numBlocks = coerce nextBlock, stateAssignment = vect }
 
 -- Could be made faster by doing what commonRefinement is doing but
 -- stopping early. This is already much faster than what is in
@@ -47,6 +51,13 @@ data Comparison
   | Incomparable
   deriving (Eq, Ord, Read, Show, Enum, Bounded)
 
+-- We put the finer partitions above
+toPreorder :: Comparison -> PartialOrdering
+toPreorder Equivalent = EQ'
+toPreorder Refinement = GT'
+toPreorder Coarsening = LT'
+toPreorder Incomparable = IC'
+
 -- See comment at isRefinementOf2
 comparePartitions :: Partition -> Partition -> Comparison
 comparePartitions p1 p2

src/Preorder.hs

@@ -0,0 +1,83 @@
+{-# language PatternSynonyms #-}
+module Preorder where
+
+{- |
+This modules includes some algorithms to deal with preorders. For our use-case
+it could be done efficiently with a single function. But this makes it a bit
+unwieldy. So I have separated it into two functions:
+1. One function takes a preorder and computes the equivalence classes.
+2. The second function takes the order into account (now only on the
+   representatives of the first function) and returns the "top" elements.
+-}
+
+import Control.Monad.Trans.Writer.Lazy (runWriter, tell)
+import Data.Bifunctor
+import Data.Foldable (find)
+import Data.Map.Strict qualified as Map
+import Data.Set qualified as Set
+
+type PartialOrdering = Maybe Ordering
+
+pattern EQ', LT', GT', IC' :: PartialOrdering
+pattern EQ' = Just EQ -- ^ Equivalent (or equal)
+pattern LT' = Just LT -- ^ Strictly less than
+pattern GT' = Just GT -- ^ Strictly greater than
+pattern IC' = Nothing -- ^ Incomparable
+
+-- | A preorder should satisfy reflexivity and transitivity. It is not assumed
+-- to be anti-symmetric.
+type Preorder x = x -> x -> PartialOrdering
+
+-- * Equivalences
+
+-- | This type captures equivalence classes by mapping each element to its
+-- representative in the same class. At the moment, the representative is
+-- absent from the map.
+-- TODO: This could be a lazy list, would that be useful?
+type EquivalenceClasses l = Map.Map l l
+
+-- | The representatives are simply returned as a list of unique elements
+type Representatives a = [a]
+
+-- | This functions takes elements of a pre-order and computes its equivalence
+-- classes. It returns one representative per class and maps each element to
+-- the representative. Number of queries to the preorder is O(n * c), where n
+-- is the number of elements and c <= n is the number of equivelance classes.
+-- The function is lazy in the second component, so if you only need to know
+-- the representatives, they can be streamed.
+equivalenceClasses :: Ord l => Preorder x -> [(l, x)] -> (EquivalenceClasses l, Representatives (l, x))
+equivalenceClasses comp ls = runWriter (go ls [] Map.empty)
+  where
+    -- end of list: return the classes
+    go [] _ classes = return classes
+    -- element x, we compare to all currently known representatives with 'find'
+    go (p@(l1, x):xs) repr classes =
+      case find (\(_, y) -> comp x y == EQ') repr of
+        -- If it is equivalent: insert in the map
+        Just (l2, _) -> go xs repr (Map.insert l1 l2 classes)
+        -- If not, add as representative, and output it
+        Nothing -> tell (pure p) >> go xs (p:repr) classes
+
+-- * Order
+
+-- | Computes the maximal elements in the preorder. Also computes the down-sets
+-- of each of those maximal elements. Note: This function expects a poset, so
+-- you should first use 'equivalenceClasses' to get representatives.
+-- Number of queries to the preorder is O(n^2). I am sure there is a better
+-- algorithm, which requires fewer comparisons. But this one is simple.
+maximalElements :: Ord l => Preorder x -> [(l, x)] -> ([(l, x)], Map.Map l [l])
+maximalElements comp ls = (maxs, downSets)
+  where
+    -- bigger first
+    ordPair p GT' = p
+    ordPair (l1, l2) LT' = (l2, l1)
+    ordPair _ _ = error "Cannot happen"
+    -- all elements in relation
+    rel = [ordPair (l1, l2) c | (l1, x1) <- ls, (l2, x2) <- ls, l1 < l2, let c = comp x1 x2, c /= IC']
+    -- elements in second component are lower
+    lows = Set.fromList . fmap snd $ rel
+    -- keep everything which is not low
+    maxs = filter (\(l, _) -> l `Set.notMember` lows) ls
+    -- keep relations from the top elements, to form a map
+    maxRel = filter (\(l, _) -> l `Set.notMember` lows) rel
+    downSets = Map.fromListWith (++) . fmap (second pure) $ maxRel