diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,30 @@ +Copyright Alexandr Sadovnikov (c) 2017 + +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + + * Redistributions in binary form must reproduce the above + copyright notice, this list of conditions and the following + disclaimer in the documentation and/or other materials provided + with the distribution. + + * Neither the name of Author name here nor the names of other + contributors may be used to endorse or promote products derived + from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -0,0 +1,111 @@ +# Math.Grads + +[![Travis](https://img.shields.io/travis/biocad/math-grads.svg)](https://travis-ci.org/biocad/math-grads) +[![hackage](https://img.shields.io/hackage/v/math-grads.svg)](https://hackage.haskell.org/package/math-grads) +[![hackage-deps](https://img.shields.io/hackage-deps/v/math-grads.svg)](https://hackage.haskell.org/package/math-grads) + +Math.Grads is library that provides graph-like data structures +and various useful algorithms for analysis of these data structures. + +Its main feature is that all of provided type classes, data structures and +functions are written in most abstract way possible to meet different demands +in functionality. + +## Data Structures + +### Graph + +Graph is a type class that upon being instantiated gives data structure +properties of graph-like object. + +### GenericGraph + +GenericGraph is a data structure that describes undirected graphs and is +parametrized by type of graph's vertices and type of graph's edges. +So it's really up to the developer what will be stored in Generic Graph's vertices +and edges. + +GenericGraph is honest instance of Graph, therefore it can be used in all functions +that require their parameters to be Graphs. + +## Algorithms + +### Ullman's subgraph isomorphism algorithm + +Math.Grads contains implementation of Ullman's subgraph isomorphism +[algorithm](https://www.cs.bgu.ac.il/~dinitz/Course/SS-12/Ullman_Algorithm.pdf). +There are several functions that one can find helpful in order to check two graphs +for isomorphism or subgraph isomorphism: + +* `isIso` checks whether two graphs are isomorphic; +* `isIsoSub` checks whether second graph has subgraph isomorphic to the first one; +* `getIso` finds matching of vertices of first graph to vertices of subgraph in second graph that +is isomorphic to the first graph; +* `getMultiIso` finds all such matchings. + +In order for these functions to work graphs that are being passed to them have to also +be instances of `GComparable` type class. + +Definition of this class is as follows: + +```haskell +class (Graph g1, Graph g2) => GComparable g1 v1 e1 g2 v2 e2 where + vComparator :: g1 v1 e1 -> g2 v2 e2 -> VComparator v1 v2 + eComparator :: g1 v1 e1 -> g2 v2 e2 -> EComparator e1 e2 + +-- | Function that checks whether two vertices are identical. +type VComparator v1 v2 = VertexIndex -> VertexIndex -> Bool + +-- | Function that checks whether two edges are identical. +type EComparator e1 e2 = GraphEdge e1 -> GraphEdge e2 -> Bool +``` + +So, basically, if two `Graph`s are `GComparable` with each other there exist +two functions that are responsible for establishing equality between vertices and edges +of such Graphs. + +Here Math.Grads gets its chance to shine, because developer isn't constrained to +what we (as developers of Math.Grads) thought would be an appropriate way for comparing +vertices and edges of your data structure. We give the developers opportunity to define +such relations for their data structures themselves. + +Maybe you want to know surroundings of two vertices in order to compare them, maybe +you don't — the choice is yours! + +### Algorithm for calculation of planar graph's coordinates + +Math.Grads provides algorithm for calculation of coordinates of planar graphs. +Its main idea is that most such graphs used in practice can be represented +as union of systems of conjugated cycles and paths that connect these systems. + +So, if you know, that your planar graph looks just like this +(for example, small molecules from chemistry perfectly fit +into the definition of graphs that can be drawn correctly by the algorithm), +you may find `getCoordsForGraph` function quite useful. + +Algorithm first draws systems of conjugated cycles, then draws paths between them, +unites systems with path and using random generator samples different conformations +of resulting graph until conformation without self-intersections (that's why graph needs +to be planar) is found. + +Once again, in order for you graph to be drawn you need to make it an instance of +special type class: + +```haskell +class Graph g => Drawable g v e where + edgeFixator :: g v e -> EdgeFixator e + edgeFixator = const $ (,) [] + +type EdgeFixator e = CoordMap -> (EdgeList e, CoordMap) +``` + +`edgeFixator` is function that given `Graph` returns other function that somehow transforms +coordinates of graph before sampling and states, which edges of graph shouldn't change their coordinates +during sampling ('fixates' them, if you will). As you can see, `edgeFixator` has default implementation, +so if you don't want such functionality, just instantiate your graph as `Drawable` without +getting into such details. + +### Miscellaneous functions on graphs + +Math.Grads also provides all other kinds of graph algorithms that you might find useful: +your depth-first searches, breadth-first searches, functions to find cycles in graphs and so on. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/math-grads.cabal b/math-grads.cabal new file mode 100644 --- /dev/null +++ b/math-grads.cabal @@ -0,0 +1,104 @@ +cabal-version: >=1.10 +name: math-grads +version: 0.1.5.1 +license: BSD3 +license-file: LICENSE +copyright: 2017 Alexandr Sadovnikov +maintainer: artemkondyukov, AlexKaneRUS, vks4git +author: Alexandr Sadovnikov +homepage: https://github.com/biocad/math-grads#readme +synopsis: Library containing graph data structures and graph algorithms +description: + Library containing graph data structures and graph algorithms. + . + Graph data structures: + . + * Graph type class; + . + * GenericGraph data structure. + . + Graph algorithms: + . + * Ullmann's subgraph isomorphism algorithm; + . + * drawing of planar graphs. +category: Math, Graph +build-type: Simple +extra-source-files: + README.md + +source-repository head + type: git + location: https://github.com/biocad/math-grads + +library + exposed-modules: + Math.Grads.Algo.Cycles + Math.Grads.Algo.Interaction + Math.Grads.Algo.Isomorphism + Math.Grads.Algo.Paths + Math.Grads.Algo.Traversals + Math.Grads.Drawing.Coords + Math.Grads.Graph + Math.Grads.GenericGraph + Math.Grads.Utils + hs-source-dirs: src + other-modules: + Math.Grads.Drawing.Internal.Coords + Math.Grads.Drawing.Internal.Cycles + Math.Grads.Drawing.Internal.CyclesPathsAlignment + Math.Grads.Drawing.Internal.Paths + Math.Grads.Drawing.Internal.Sampling + Math.Grads.Drawing.Internal.Utils + Math.Grads.Angem + Math.Grads.Angem.Internal.VectorOperations + Math.Grads.Angem.Internal.MatrixOperations + default-language: Haskell2010 + build-depends: + base >=4.7 && <5, + aeson <1.5, + array <0.6, + containers <0.7, + linear <1.21, + matrix <0.4, + mtl <2.3, + random <1.2, + vector <0.13 + +test-suite Coords-test + type: exitcode-stdio-1.0 + main-is: Coords.hs + hs-source-dirs: test + default-language: Haskell2010 + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base <4.13, + containers <0.7, + hspec <2.7, + math-grads -any, + random <1.2 + +test-suite Graph-test + type: exitcode-stdio-1.0 + main-is: Graph.hs + hs-source-dirs: test + default-language: Haskell2010 + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base <4.13, + containers <0.7, + hspec <2.7, + math-grads -any + +test-suite Isomorphism-test + type: exitcode-stdio-1.0 + main-is: Isomorphism.hs + hs-source-dirs: test + default-language: Haskell2010 + ghc-options: -threaded -rtsopts -with-rtsopts=-N + build-depends: + base <4.13, + array <0.6, + containers <0.7, + hspec <2.7, + math-grads -any diff --git a/src/Math/Grads/Algo/Cycles.hs b/src/Math/Grads/Algo/Cycles.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Algo/Cycles.hs @@ -0,0 +1,155 @@ +-- | Module that provides functions for analysis of graph's cycles. +-- +module Math.Grads.Algo.Cycles + ( findCycles + , findLocalCycles + , getCyclic + , isEdgeInCycle + ) where + +import Control.Monad.State (State, runState) +import Control.Monad.State.Class (get, modify) +import Data.List (partition, sort, union, (\\)) +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M (keys, (!)) +import Data.Set (Set) +import qualified Data.Set as S (empty, fromList, insert, + member) +import Math.Grads.Algo.Interaction (edgeListToMap, getEnds, + getIndices, getOtherEnd, + getVertexIncident, haveSharedEdge, + matchEdges) +import Math.Grads.Algo.Paths (dfsAllPaths) +import Math.Grads.GenericGraph (GenericGraph, safeIdx) +import Math.Grads.Graph (EdgeList, GraphEdge, vCount) + +-- | Takes 'EdgeList' and finds non-redundant set of conjugated simple cycles. +-- Cycles sharing in common one edge are considered to be one cycle. +-- BondList must obey rule (b, e, _) b < e. +-- +findCycles :: Ord e => EdgeList e -> [EdgeList e] +findCycles bonds = sort <$> conjRings redundantCycles + where + redundantCycles = findCyclesR bonds + + findCyclesR :: Ord e => EdgeList e -> [EdgeList e] + findCyclesR bs = let (result, taken) = stateCycles bs in + if sort taken == sort bs then result + else result ++ findCyclesR (bs \\ taken) + + stateCycles :: Ord e => EdgeList e -> ([EdgeList e], EdgeList e) + stateCycles bs = runState (cyclesHelper bs [] (minimum (getIndices bs))) [] + +conjRings :: Ord e => [EdgeList e] -> [EdgeList e] +conjRings (b : bs) = + let + (shd, rest) = partition (haveSharedEdge b) bs + in + case shd of + [] -> b : conjRings rest + _ -> conjRings $ foldr union b shd : rest +conjRings b = b + +takeCycle :: EdgeList e -> GraphEdge e -> EdgeList e +takeCycle [] _ = [] +takeCycle bl@((aPop, bPop, _) : _) bn@(aNow, bNow, _) = bn : takeWhile cond bl ++ take 1 (dropWhile cond bl) + where + theB | bPop == aNow = bNow + | bPop == bNow = aNow + | aPop == bNow = aNow + | otherwise = bNow + cond :: GraphEdge e -> Bool + cond (a', b', _) = theB /= a' && theB /= b' + +cyclesHelper :: Eq e => EdgeList e -> EdgeList e -> Int -> State (EdgeList e) [EdgeList e] +cyclesHelper bs trc n = do + curSt <- get + let adjBonds = filter (`notElem` curSt) $ getVertexIncident bs n + + let visited = concatMap getEnds curSt + let curBondClosures = filter (\b -> getOtherEnd b n `elem` visited) adjBonds + let furtherBonds = filter (`notElem` curBondClosures) adjBonds + + let procBnd bnd = cyclesHelper bs (bnd : trc) (getOtherEnd bnd n) + restBondClosures <- mapM (\b -> modify (b:) >>= const (procBnd b)) furtherBonds + + return $ (takeCycle trc <$> curBondClosures) ++ concat restBondClosures + +-- | Checks that edge with given index in 'EdgeList' is contained in any cycle. +-- +isEdgeInCycle :: Ord e => EdgeList e -> Int -> Bool +isEdgeInCycle bs n = any ((bs !! n) `elem`) $ findCycles bs + +-- | Finds all cycles of minimal length contained in system of conjugated cycles. +-- +findLocalCycles :: Eq e => EdgeList e -> [EdgeList e] +findLocalCycles bonds = if null cycles then [] + else helperFilter (tail res) [head res] + where + -- TODO: We need to remove this filter. + cycles = filter (\x -> length x < 21) (findLocalCycles' bonds) + res = filter (`filterBigCycles` cycles) cycles + +findLocalCycles' :: Eq e => EdgeList e -> [EdgeList e] +findLocalCycles' bonds = concatMap (\(a, b, _) -> dfsAllPaths bonds a b) cycleBonds + where + stBonds = dfsSt bonds + cycleBonds = bonds \\ stBonds + + dfsSt :: EdgeList e -> EdgeList e + dfsSt bs = matchEdges bs bondsInd + where + graph = edgeListToMap bs + bondsInd = dfsSt' graph (M.keys graph) [] [] + + dfsSt' :: Map Int [Int] -> [Int] -> [Int] -> [(Int, Int)] -> [(Int, Int)] + dfsSt' _ [] _ bs = bs + dfsSt' graph (current : toVisit) visited bs | current `elem` visited = dfsSt' graph toVisit visited bs + | otherwise = dfsSt' graph toVisitModified (current:visited) visitedBonds + where + visitedBonds = bs ++ if not (null visited) then [found | snd found /= -1] else [] + found = findRib graph visited current + + toVisitModified = (graph M.! current) ++ toVisit + + findRib :: Map Int [Int] -> [Int] -> Int -> (Int, Int) + findRib graph visited current = (current, if not (null found) then head found else -1) + where + found = filter (`elem` visited) (graph M.! current) + +filterBigCycles :: Eq e => EdgeList e -> [EdgeList e] -> Bool +filterBigCycles currentCycle cycles = not (foldl (\x y -> x || currentCycle /= y && length currentCycle > length y && length (filter (`elem` currentCycle) y) > 1) False cycles) + +helperFilter :: Eq e => [EdgeList e] -> [EdgeList e] -> [EdgeList e] +helperFilter [] ready = ready +helperFilter (x:xs) ready = if exists x ready then helperFilter xs ready else helperFilter xs (x:ready) + where + exists a1 = any (\x' -> length a1 == length x' && all (\(a, b, t) -> (a, b, t) `elem` x' || (b, a, t) `elem` x') a1) + +-- | Checks whether or not given vertex belongs to any cycle. +-- +isCyclic :: GenericGraph v e -> Int -> Int -> (Bool, Set Int) -> Int -> (Bool, Set Int) +isCyclic graph target previous (result, visited) current | result = (result, visited) + | (previous /= (-1)) && (current == target) = (True, visited) + | current `S.member` visited = (result, visited) + | otherwise = foldl foldFunc (result, updatedVis) next + where + next :: [Int] + next = filter (/= previous) $ graph `safeIdx` current + + updatedVis :: Set Int + updatedVis = current `S.insert` visited + + foldFunc :: (Bool, Set Int) -> Int -> (Bool, Set Int) + foldFunc = isCyclic graph target current + +-- | Returns the set of all vertices which belong to any cycle. +-- +getCyclic :: GenericGraph v e -> Set Int +getCyclic graph = S.fromList . map fst . filter snd $ zip indices cyclic + where + indices :: [Int] + indices = [0 .. vCount graph - 1] + + cyclic :: [Bool] + cyclic = map (\ix -> fst $ isCyclic graph ix (-1) (False, S.empty) ix) indices diff --git a/src/Math/Grads/Algo/Interaction.hs b/src/Math/Grads/Algo/Interaction.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Algo/Interaction.hs @@ -0,0 +1,149 @@ +-- | Module that provides various functions for interaction with 'GraphEdge's, +-- 'EdgeList's and vertices themselves. +-- +module Math.Grads.Algo.Interaction + ( + -- * Vertex Functions + -- + areAdjacent + , getEnds + , getOtherEnd + , getSharedVertex + , getVertexAdjacent + , getVertexIncident + , getVertexIncidentIdx + , haveSharedVertex + , isIncident + , (~=) + , (/~=) + -- * Edge Functions + -- + , matchEdges + , getEdgeIncident + -- * EdgeList Functions + -- + , doubleEdgeList + , edgeListToMap + , haveSharedEdge + , sortBondList + , getIndices + ) where + +import Data.List (find, findIndices, intersect, sortOn) +import Data.Map (Map) +import qualified Data.Map as M +import Data.Maybe (fromJust, isJust) + +import Math.Grads.Graph (EdgeList, GraphEdge, edgeType) +import Math.Grads.Utils (nub) + +-- | Equality operator for 'GraphEdge's. +-- +(~=) :: GraphEdge e1 -> GraphEdge e2 -> Bool +(b, e, _) ~= (b', e', _) = (b == b' && e == e') || (b == e' && e == b') + +-- | Inequality operator for 'GraphEdge's. +-- +(/~=) :: GraphEdge e1 -> GraphEdge e2 -> Bool +b1 /~= b2 = not $ b1 ~= b2 + +-- | Checks that vertex with given index is incident to 'GraphEdge'. +-- +isIncident :: GraphEdge e -> Int -> Bool +isIncident (b, e, _) n = b == n || e == n + +-- | Find all edges in given 'EdgeList' that are incident to vertex with given index. +-- +getVertexIncident :: EdgeList e -> Int -> EdgeList e +getVertexIncident bs n = filter (`isIncident` n) bs + +-- | Returns indices of edges in 'EdgeList' that are incident to vertex with given index. +-- +getVertexIncidentIdx :: EdgeList e -> Int -> [Int] +getVertexIncidentIdx bs n = findIndices (`isIncident` n) bs + +-- | Returns index of vertex incident to given 'GraphEdge' and different from passed index. +-- +getOtherEnd :: GraphEdge e -> Int -> Int +getOtherEnd (b, e, _) n | b == n = e + | e == n = b + | otherwise = error "There is no such index in edge." + +-- | Finds in given 'EdgeList' all indices of vertices adjacent to given vertex. +-- +getVertexAdjacent :: EdgeList e -> Int -> [Int] +getVertexAdjacent bs n = (`getOtherEnd` n) <$> getVertexIncident bs n + +-- | Checks whether two vertices with given indices are adjacent in given 'EdgeList'. +-- +areAdjacent :: EdgeList e -> Int -> Int -> Bool +areAdjacent bs n n' = n' `elem` getVertexAdjacent bs n + +-- | Retrieves indices of vertices that are being connected by given 'GraphEdge'. +-- +getEnds :: GraphEdge e -> [Int] +getEnds (b, e, _) = [b, e] + +-- | Checks that two edges have common vertex. +-- +haveSharedVertex :: GraphEdge e1 -> GraphEdge e2 -> Bool +haveSharedVertex b1 b2 = isJust $ getSharedVertex b1 b2 + +-- | Gets shared common vertex of two edges. If edges don't have common vertex, +-- returns Nothing. +-- +getSharedVertex :: GraphEdge e1 -> GraphEdge e2 -> Maybe Int +getSharedVertex b1 b2 | null is = Nothing + | length is == 2 = Nothing + | otherwise = Just $ head is + where + is = getEnds b1 `intersect` getEnds b2 + +-- | Find edges in 'EdgeList' which ordered pairs of indices, that they are connecting, +-- are present in passed list of ordered pairs. +-- +matchEdges :: EdgeList e -> [(Int, Int)] -> EdgeList e +matchEdges bonds = fmap (\(a, b) -> fromJust (find (~= (a, b, undefined)) bonds)) + +-- | Find all edges that are incident to edge in 'EdgeList' with given index. +-- +getEdgeIncident :: Ord e => EdgeList e -> Int -> EdgeList e +getEdgeIncident bs n | n >= length bs = [] + | otherwise = filter (/= (beg, end, typ)) $ getVertexIncident bs beg ++ getVertexIncident bs end + where + (beg, end, typ) = bs !! n + +-- | For every edge in 'EdgeList' add to that list an edge in opposite direction. +-- +doubleEdgeList :: EdgeList e -> EdgeList e +doubleEdgeList = concatMap (\(a, b, t) -> [(a, b, t), (b, a, t)]) + +-- | Transforms 'EdgeList' into 'Map' that corresponds to adjacency list of undirected +-- graph induced by these edges. +-- +edgeListToMap :: EdgeList e -> Map Int [Int] +edgeListToMap bonds' = M.fromList (fmap (toVertex bonds) (getIndices bonds)) + where + bonds = doubleEdgeList bonds' + + toVertex :: EdgeList e -> Int -> (Int, [Int]) + toVertex bs i = (i, concatMap (\(a, b, _) -> [a | b == i]) bs) + +-- | Checks that two 'EdgeList's have common edge. +-- +haveSharedEdge :: Eq e => EdgeList e -> EdgeList e -> Bool +haveSharedEdge b1 b2 = or $ fmap (`elem` b2) b1 + +-- | Sorting for 'EdgeList', that sorts edges on their type, then on index of their +-- to (right) vertex, then on index of their from (left) vertex. +-- +sortBondList :: Ord e => EdgeList e -> EdgeList e +sortBondList = sortOn left . sortOn right . sortOn edgeType + where + left (a, _, _) = a + right (_, b, _) = b + +-- | Gets all vertices from 'EdgeList'. +-- +getIndices :: EdgeList e -> [Int] +getIndices = nub . concatMap getEnds diff --git a/src/Math/Grads/Algo/Isomorphism.hs b/src/Math/Grads/Algo/Isomorphism.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Algo/Isomorphism.hs @@ -0,0 +1,301 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Module that provides functions for identifying graph/subgraph isomorphism. +-- +module Math.Grads.Algo.Isomorphism + ( EComparator, VComparator + , GComparable (..) + , VertexIndex + , getIso + , getMultiIso + , isIso + , isIsoSub + ) where + +import Control.Arrow (second, (&&&), (***)) +import qualified Data.Array as A +import Data.List (delete, sortOn) +import Data.Map (Map) +import qualified Data.Map.Strict as M +import Data.Matrix (Matrix (..), getElem, getRow, mapRow, + matrix, multStd, ncols, nrows, + setElem, transpose) +import Data.Maybe (isJust, listToMaybe) +import Data.Tuple (swap) +import qualified Data.Vector as V +import Math.Grads.GenericGraph (GenericGraph (..)) +import Math.Grads.Graph (Graph, GraphEdge, changeIndsEdge, + fromList, incidentIdx, toList, vCount, + (!.)) +import Math.Grads.Utils (nub) + +type GenericGraphIso v e = GenericGraph Int e + +-- | Type alias for 'Int'. +-- +type VertexIndex = Int + +-- | Function that checks whether two vertices are identical. +-- Due to properties related to index of vertex, +-- like number of neighbors, we consider vertex indices instead of vertices. +-- +type VComparator v1 v2 = VertexIndex -> VertexIndex -> Bool + +-- | Function that checks whether two edges are identical. +-- Due to properties related to index of vertex, +-- like belonging to a cycle, we consider GraphEdge (Int, Int, e) instead of e. +-- +type EComparator e1 e2 = GraphEdge e1 -> GraphEdge e2 -> Bool + +-- | Type class for graphs that could be checked for isomorphism. +-- +class (Graph g1, Graph g2) => GComparable g1 v1 e1 g2 v2 e2 where + vComparator :: g1 v1 e1 -> g2 v2 e2 -> VComparator v1 v2 + eComparator :: g1 v1 e1 -> g2 v2 e2 -> EComparator e1 e2 + +-- | Checks whether two graphs are isomorphic. +-- +isIso :: (Ord v1, Ord v2, GComparable GenericGraph v1 e1 GenericGraph v2 e2, Eq e1, Eq e2) + => GenericGraph v1 e1 + -> GenericGraph v2 e2 + -> Bool +isIso queryGraph targetGraph = res + where + (v1, e1) = toList queryGraph + (v2, e2) = toList targetGraph + isoSub = isIsoSub queryGraph targetGraph + + res = length v1 == length v2 && length e1 == length e2 && isoSub + +-- | Check for queryGraph \( \subseteq \) targetGraph. +-- +isIsoSub :: (Ord v1, Ord v2, GComparable GenericGraph v1 e1 GenericGraph v2 e2, Eq e1, Eq e2) + => GenericGraph v1 e1 -- ^ queryGraph + -> GenericGraph v2 e2 -- ^ targetGraph + -> Bool +isIsoSub queryGraph targetGraph = isJust $ getIso queryGraph targetGraph + +-- | Get one vertices matching (if exists) from queryGraph to targetGraph. +-- +getIso :: (Ord v1, Ord v2, GComparable GenericGraph v1 e1 GenericGraph v2 e2, Eq e1, Eq e2) + => GenericGraph v1 e1 -- ^ queryGraph + -> GenericGraph v2 e2 -- ^ targetGraph + -> Maybe (Map Int Int) +getIso queryGraph targetGraph = listToMaybe $ getMultiIso queryGraph targetGraph + +-- | Get all possible vertices matchings from queryGraph to targetGraph. +-- +getMultiIso :: (Ord v1, Ord v2, GComparable GenericGraph v1 e1 GenericGraph v2 e2, Eq e1, Eq e2) + => GenericGraph v1 e1 -- ^ queryGraph + -> GenericGraph v2 e2 -- ^ targetGraph + -> [Map Int Int] +getMultiIso queryGraph' targetGraph' = matches + where + ((queryGraph, queryGraphWI), fromIsoToOldQ) = second inverseMap $ graphToGraphIso queryGraph' + ((targetGraph, targetGraphWI), fromIsoToOldT) = second inverseMap $ graphToGraphIso targetGraph' + + vComp = vComparator queryGraphWI targetGraphWI + eComp = eComparator queryGraphWI targetGraphWI + + isos = isoGraph vComp eComp queryGraph targetGraph + matches = fmap (\x -> getMatchMap x fromIsoToOldQ fromIsoToOldT) isos + +inverseMap :: Map Int Int -> Map Int Int +inverseMap = M.fromList . (swap <$>) . M.toList + +getMatchMap :: Matrix Int -> Map Int Int -> Map Int Int -> Map Int Int +getMatchMap isoMatrix fromIsoToOldQ fromIsoToOldT = res + where + forMap = fmap (getMatchRow isoMatrix) [0 .. nrows isoMatrix - 1] + res = M.fromList (fmap ((fromIsoToOldQ M.!) *** (fromIsoToOldT M.!)) forMap) + +getMatchRow :: Matrix Int -> Int -> (Int, Int) +getMatchRow isoMatrix ind = (ind, helper 0) + where + row = getRow (ind + 1) isoMatrix + + helper :: Int -> Int + helper counter = if row V.! counter == 1 then counter + else helper (counter + 1) + +isoGraph :: (Eq e1, Eq e2) => VComparator v1 v2 + -> EComparator e1 e2 + -> GenericGraphIso v1 e1 + -> GenericGraphIso v2 e2 + -> [Matrix Int] +isoGraph vComp eComp queryGraph targetGraph = res + where + queryGraphEdges = (fst <$>) <$> gAdjacency queryGraph + sizeOfQueryGraph = vCount queryGraph + pMatrix = matrix sizeOfQueryGraph sizeOfQueryGraph (\(i, j) -> if i - 1 `elem` queryGraphEdges A.! (j - 1) then 1 else 0) + + targetGraphEdges = (fst <$>) <$> gAdjacency targetGraph + sizeOfTargetGraph = vCount targetGraph + gMatrix = matrix sizeOfTargetGraph sizeOfTargetGraph (\(i, j) -> if i - 1 `elem` targetGraphEdges A.! (j - 1) then 1 else 0) + + mMatrix = matrix sizeOfQueryGraph sizeOfTargetGraph (\(i, j) -> if fits vComp eComp queryGraph targetGraph i j then 1 else 0) + + currentRow = 0 + unusedColumns = [1 .. ncols mMatrix] + + res = recurse eComp queryGraph targetGraph unusedColumns currentRow gMatrix pMatrix mMatrix + +fits :: (Eq e1, Eq e2) => VComparator v1 v2 + -> EComparator e1 e2 + -> GenericGraphIso v1 e1 + -> GenericGraphIso v2 e2 + -> Int + -> Int + -> Bool +fits vComp eComp queryGraph targetGraph i j = res + where + (vertex, edges) = (gIndex queryGraph A.! (i - 1), incidentIdx queryGraph $ i - 1) + (vertex', edges') = (gIndex targetGraph A.! (j - 1), incidentIdx targetGraph $ j - 1) + res = length edges <= length edges' && canBeSubset eComp edges edges' && vertex `vComp` vertex' + +canBeSubset :: forall e1 e2. EComparator e1 e2 -> [GraphEdge e1] -> [GraphEdge e2] -> Bool +canBeSubset eComp query target = uniqueSeq maps + where + bondsInd = zip [0 ..] target + maps = findMatches <$> query + + findMatches :: GraphEdge e1 -> [Int] + findMatches thisEdge = fst <$> filter (\(_, otherEdge) -> eComp thisEdge otherEdge) bondsInd + +uniqueSeq :: [[Int]] -> Bool +uniqueSeq maps = res + where + seqs = sequence maps + + res = any (\x -> length x == length (nub x)) seqs + +-- | Converts input graph into graph in which vertices with most amount of edges have lowest indices. +-- +graphToGraphIso :: (Ord v) => GenericGraph v e + -> ((GenericGraphIso v e, GenericGraph v e), M.Map Int Int) +graphToGraphIso graph = res + where + (vertices, edges) = toList graph + vArr = gIndex graph + + indsWithNCount = fmap (id &&& (length . (graph !.))) [0 .. length vertices - 1] + sortedInds = fst <$> sortOn (\x -> - (snd x)) indsWithNCount + changesMap = M.fromList (zip sortedInds [0 ..]) + + sortedV = fmap (vArr A.!) sortedInds + changedEdges = fmap (changeIndsEdge (changesMap M.!)) edges + + forGraphWI = (sortedV, changedEdges) + forGraph = ([0 .. length sortedV - 1], changedEdges) + + res = ((fromList forGraph, fromList forGraphWI), changesMap) + +-- | Ullman's subgraph isomorphism algorithm itself. +-- +recurse :: (Eq e1, Eq e2) => EComparator e1 e2 + -> GenericGraphIso v1 e1 + -> GenericGraphIso v2 e2 + -> [Int] + -> Int + -> Matrix Int + -> Matrix Int + -> Matrix Int + -> [Matrix Int] +recurse eComp queryGraph targetGraph unusedColumns currentRow gMatrix pMatrix mMatrix = res + where + prunedM = prune eComp queryGraph targetGraph mMatrix currentRow + + recs = concatMap pruneNext unusedColumns + + res | hasEmptyRow mMatrix = [] + | currentRow == nrows mMatrix && isIsomorphism gMatrix pMatrix mMatrix = [mMatrix] + | not (hasEmptyRow prunedM) = recs + | otherwise = [] + + pruneNext :: Int -> [Matrix Int] + pruneNext x = recurse eComp queryGraph targetGraph newColumns newRow gMatrix pMatrix changedMatrix + where + newColumns = delete x unusedColumns + newRow = currentRow + 1 + changedMatrix = changeRow prunedM newRow x + +prune :: (Eq e1, Eq e2) => EComparator e1 e2 + -> GenericGraphIso v1 e1 + -> GenericGraphIso v2 e2 + -> Matrix Int + -> Int + -> Matrix Int +prune eComp queryGraph targetGraph mMatrix currentRow | null indicesToChange = mMatrix + | hasEmptyRow mMatrix = mMatrix + | otherwise = res + where + numberOfMRows = nrows mMatrix + numberOfMColumns = ncols mMatrix + pairsOfindices = [(i, j) | i <- [1.. numberOfMRows], j <- [1.. numberOfMColumns], getElem i j mMatrix == 1] + + suitPair :: Int -> Int -> Bool + suitPair = hasSuitableNeighbors eComp queryGraph targetGraph mMatrix + + indicesToChange = filter (not . uncurry suitPair) pairsOfindices + changedMMatrix = foldl (flip (setElem 0)) mMatrix indicesToChange + + res = prune eComp queryGraph targetGraph changedMMatrix currentRow + +-- | Returns True if we can map all neighbors of query vertex to neighbors of target vertex in mMatrix. +-- +hasSuitableNeighbors :: forall v1 v2 e1 e2. (Eq e1, Eq e2) => EComparator e1 e2 + -> GenericGraphIso v1 e1 + -> GenericGraphIso v2 e2 + -> Matrix Int + -> Int + -> Int + -> Bool +hasSuitableNeighbors eComp queryGraph targetGraph mMatrix query target = doesSatisfy + where + iQ = query - 1 + iT = target - 1 + + neighborsOfQ = (\(i, e) -> (iQ, i, e)) <$> queryGraph !. iQ + neighborsOfT = (\(i, e) -> (iT, i, e)) <$> targetGraph !. iT + + hasProperNeighbor :: GraphEdge e1 -> Bool + hasProperNeighbor edge = any (\edge' -> getProperElem edge edge' == 1 && eComp edge edge') neighborsOfT + + getProperElem :: GraphEdge e1 -> GraphEdge e2 -> Int + getProperElem (_, b, _) (_, b', _) = getElem (b + 1) (b' + 1) mMatrix + + doesSatisfy = all hasProperNeighbor neighborsOfQ + +-- | Checks whether mMatrix encodes an isomorphism between pMatrix and gMatrix. +-- +isIsomorphism :: Matrix Int -- ^ gMatrix + -> Matrix Int -- ^ pMatrix + -> Matrix Int -- ^ mMatrix + -> Bool +isIsomorphism gMatrix pMatrix mMatrix = leqMatrices pMatrix check + where + check = multStd mMatrix (transpose (multStd mMatrix gMatrix)) + +-- | Componentwise "less or equal" operation for matrices. +-- +leqMatrices :: Matrix Int -> Matrix Int -> Bool +leqMatrices matrixA matrixB = nrows matrixA * ncols matrixA <= nrows matrixB * ncols matrixB && helper elems + where + numOfRows = nrows matrixA + numOfColumns = ncols matrixB + elems = [(i, j) | i <- [1..numOfRows], j <- [1..numOfColumns]] + helper = foldr (\x -> (&&) (uncurry getElem x matrixA <= uncurry getElem x matrixB)) True + +-- | Replace all elements in row with 0 apart from chosen one. +-- +changeRow :: Matrix Int -> Int -> Int -> Matrix Int +changeRow mMatrix row column = mapRow helper row mMatrix + where helper column' a = if column' /= column then 0 else a + +hasEmptyRow :: Matrix Int -> Bool +hasEmptyRow prunedMatrix = cond + where + numberOfRows = nrows prunedMatrix + cond = any (\x -> all (== 0) (getRow x prunedMatrix)) [1 .. numberOfRows] diff --git a/src/Math/Grads/Algo/Paths.hs b/src/Math/Grads/Algo/Paths.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Algo/Paths.hs @@ -0,0 +1,118 @@ +-- | Module that provides functions for different kinds of path-finding in graph. +-- +module Math.Grads.Algo.Paths + ( allPathsInGraph + , allPathsFromVertex + , dfsAllPaths + , dfsSearch + , findBeginnings + ) where + +import Control.Monad (forM_) +import Control.Monad.State (State, execState) +import Control.Monad.State.Class (get, modify) +import Data.Map (Map, keys, (!)) +import Data.Maybe (fromMaybe, isJust) +import Math.Grads.Algo.Interaction (edgeListToMap, getVertexAdjacent, + matchEdges, sortBondList) +import Math.Grads.GenericGraph (GenericGraph) +import Math.Grads.Graph (EdgeList, Graph (..)) +import Math.Grads.Utils (nub, subsets, uniter) + +-- | Finds all vertices in 'EdgeList' that have only one neighbour. +-- +findBeginnings :: EdgeList e -> [Int] +findBeginnings edges = fmap fst (filter ((== 1) . snd) counters) + where + graph = edgeListToMap edges + counters = zip (keys graph) (fmap (length . (graph !)) (keys graph)) + +-- | Calculates all branched paths in graph up to the given length. +-- +allPathsInGraph :: Ord e => GenericGraph v e -> Int -> [EdgeList e] +allPathsInGraph graph lengthOfPath = helper graph vertexInds [] + where + vertexInds = [0 .. (vCount graph - 1)] + + helper :: Ord e => GenericGraph v e -> [Int] -> [Int] -> [EdgeList e] + helper _ [] _ = [] + helper gr (x : xs) forbidden = allPathsFromVertex gr x lengthOfPath forbidden ++ helper gr xs (x : forbidden) + +-- | Calculates all branched paths up to the given length from given vertex in graph +-- considering indices of vertices that shouldn't be visited during path-finding. +-- +allPathsFromVertex :: Ord e => GenericGraph v e -> Int -> Int -> [Int] -> [EdgeList e] +allPathsFromVertex graph vertex lengthOfPath forbidden = nub filtered + where + res' = execState (allPathsFromVertexSt graph [vertex] lengthOfPath forbidden []) [] + filtered = sortBondList <$> filter (not . null) res' + +allPathsFromVertexSt :: Ord e => GenericGraph v e -> [Int] -> Int -> [Int] -> EdgeList e -> State [EdgeList e] [EdgeList e] +allPathsFromVertexSt graph vertices lenOfPath forbidden res = if lenOfPath < 0 then get + else + do + modify (res :) + + let edgesNeigh = nub (filter (`notElem` res) (concatMap (incidentIdx graph) vertices)) + let allowedEdgesNeigh = filter (\(a, b, _) -> a `notElem` forbidden && b `notElem` forbidden) edgesNeigh + let edgeSets = filter ((\x -> x > 0 && x <= lenOfPath) . length) (subsets allowedEdgesNeigh) + if lenOfPath == 0 || not (null allowedEdgesNeigh) then + do + forM_ edgeSets (\set -> do + let newNeighbors = concatMap (getVertexAdjacent set) vertices + let newLength = lenOfPath - length set + let newRes = res ++ set + modify (execState (allPathsFromVertexSt graph newNeighbors newLength forbidden newRes) [] ++)) + get + else get + +-- | Finds path between two vertices in graph represented as 'EdgeList'. +-- Graph shouldn't have any cycles. Hmmm, what's the difference between this function +-- and DFS or BFS?.. +-- +dfsSearch :: EdgeList e -> Int -> Int -> Maybe (EdgeList e, [Int]) +dfsSearch edges start finish = if cond then Just (matchEdges edges edgesInd, x) + else Nothing + where + graph = edgeListToMap edges + x = fromMaybe [] $ helperDfs graph (-1) finish [start] + + edgesInd = uniter x + inds = concatMap (\(x', y, _) -> [x', y]) edges + + cond = start `elem` inds && finish `elem` inds + +helperDfs :: Map Int [Int] -> Int -> Int -> [Int] -> Maybe [Int] +helperDfs graph prev finish path | current /= prev && current /= finish = if not (null (==?)) then head (==?) else Nothing + | current == finish = Just path + | otherwise = Nothing + where + current = head path + children = filter (/= prev) (graph ! current) + (==?) = filter isJust (map (\x -> helperDfs graph current finish (x : path)) children) + +-- | Finds all paths between vertices with given indices in 'EdgeList'. +-- +dfsAllPaths :: EdgeList e -> Int -> Int -> [EdgeList e] +dfsAllPaths edges start finish = fmap (matchEdges edges) edgesInd + where + graph = edgeListToMap edges + paths = execState (statePaths graph finish [start]) [] + + filteredPaths = filter ((> 2) . length) paths + edgesInd = fmap helper filteredPaths + + helper :: [Int] -> [(Int, Int)] + helper l = if (start, finish) `elem` united then united + else (start, finish) : united + where + united = uniter l + +statePaths :: Map Int [Int] -> Int -> [Int] -> State [[Int]] [[Int]] +statePaths graph finish path = if head path `elem` tail path then get else (do + let current = head path + if current == finish then do {modify ([path] ++); get} else + do + let children = filter (`notElem` path) (graph ! current) + forM_ children (\child -> modify (execState (statePaths graph finish (child : path)) [] ++)) + get) diff --git a/src/Math/Grads/Algo/Traversals.hs b/src/Math/Grads/Algo/Traversals.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Algo/Traversals.hs @@ -0,0 +1,98 @@ +-- | Module providing various kinds of graph traversals and their modifications. +-- +module Math.Grads.Algo.Traversals + ( BFSState + , bfsState + , dfsCycle + , dfs + , getComps + , getCompsIndices + ) where + +import Control.Arrow ((&&&)) +import Control.Monad.State (State, execState) +import Control.Monad.State.Class (get, put) +import qualified Data.Array as A +import Data.List (findIndex) +import Data.Map (Map, keys, (!)) +import Data.Maybe (fromJust) + +import Math.Grads.Algo.Interaction (edgeListToMap, getIndices, + getOtherEnd, matchEdges, (~=)) +import Math.Grads.GenericGraph (GenericGraph, gIndex, subgraph) +import Math.Grads.Graph (EdgeList, Graph (..)) +import Math.Grads.Utils (nub) + +-- | Classic dfs. +-- +dfs :: EdgeList e -> Int -> EdgeList e +dfs bonds ind = if ind `elem` keys graphMap then matchEdges bonds bondsInd else [] + where + graphMap = edgeListToMap bonds + bondsInd = dfs' graphMap [ind] [] [] + +dfs' :: Map Int [Int] -> [Int] -> [Int] -> [(Int, Int)] -> [(Int, Int)] +dfs' _ [] _ bonds = bonds +dfs' gr (cur : rest) vis bs | cur `elem` vis = dfs' gr rest vis bs + | otherwise = dfs' gr (gr ! cur ++ rest) (cur : vis) visitedBonds + where + visitedBonds = concatMap helper (gr ! cur) ++ bs + + helper :: Int -> [(Int, Int)] + helper sec = [(cur, sec) | notElem (cur, sec) bs && notElem (sec, cur) bs] + +-- | Get connected components of graph. +-- Note that indexation will be CHANGED. +-- +getComps :: Ord v => GenericGraph v e -> [GenericGraph v e] +getComps graph = res + where + (_, edges) = toList graph + comps = getComps' edges [0..length (gIndex graph) - 1] [] [] + res = fmap (subgraph graph) comps + +-- | Get indices of vertices that belong to different connected components. +-- +getCompsIndices :: Ord v => GenericGraph v e -> [[Int]] +getCompsIndices graph = getComps' (snd $ toList graph) [0..length (gIndex graph) - 1] [] [] + +getComps' :: EdgeList e -> [Int] -> [Int] -> [[Int]] -> [[Int]] +getComps' _ [] _ res = res +getComps' edges (x : xs) taken res = if x `elem` taken then getComps' edges xs taken res + else getComps' edges xs (taken ++ newComp) (newComp : res) + where + newComp = nub (x : getIndices (dfs edges x)) + +-- | Dfs a simple cycle. +-- +dfsCycle :: A.Array Int [Int] -> [Int] -> [Int] -> [Int] +dfsCycle _ [] visited = visited +dfsCycle graph (current:toVisit) visited | current `elem` visited = dfsCycle graph toVisit visited + | otherwise = dfsCycle graph ((graph A.! current) ++ toVisit) (current:visited) + +-- | List of (level, (edgeIdx, vertexIdx)). +-- +type BFSState = [(Int, (Int, Int))] + +-- | Bfs algorithm that takes graph, its 'EdgeList', 'BFSState' corresponding to +-- already visited vertices and 'BFSState' that corresponds to starting point +-- of traversal and returns 'BFSState' as a result. +-- +bfsState :: (Ord v, Eq e, Show v, Show e, Graph g) => g v e -> EdgeList e -> BFSState -> BFSState -> BFSState +bfsState graph bonds ign start = fst $ execState (bfsState' graph bonds) (ign, start) + +-- | Traverses graph from a given starting point (queue) in Breadth-first search manner. +-- +bfsState' :: (Ord v, Eq e, Show v, Show e, Graph g) => g v e -> EdgeList e -> State (BFSState, BFSState) () +bfsState' gr bonds = do + (visited, queue) <- get + let (visitedL, (visitedB, visitedV)) = (fst &&& unzip . snd) $ unzip visited + case queue of + ((curLevel, (curBnd, curNum)) : rest) -> do + let curInc = (fromJust . (\x -> (~= x) `findIndex` bonds)) <$> gr `incidentIdx` curNum + let nextBonds = nub $ filter ((`notElem` visitedV) . (`getOtherEnd` curNum) . (bonds !!)) curInc + let nextLevel = ((`getOtherEnd` curNum) . (bonds !!)) <$> nextBonds + let nextVisited = zip (curLevel : visitedL) $ zip (curBnd : visitedB) (curNum : visitedV) + put (nextVisited, rest ++ zip (repeat $ curLevel +1) (zip nextBonds nextLevel)) + bfsState' gr bonds + _ -> return () diff --git a/src/Math/Grads/Angem.hs b/src/Math/Grads/Angem.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Angem.hs @@ -0,0 +1,15 @@ +-- | Module providing miscellaneous functions for working with +-- coordinates, vectors and matrices. +-- +module Math.Grads.Angem + ( alignmentFunc + , areIntersected + , eqV2 + , rotation2D + , reflectPoint + ) where + +import Math.Grads.Angem.Internal.MatrixOperations (alignmentFunc, + rotation2D) +import Math.Grads.Angem.Internal.VectorOperations (areIntersected, + eqV2, reflectPoint) diff --git a/src/Math/Grads/Angem/Internal/MatrixOperations.hs b/src/Math/Grads/Angem/Internal/MatrixOperations.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Angem/Internal/MatrixOperations.hs @@ -0,0 +1,89 @@ +-- | Functions for working with coordinates' alignment and matrix rotations. +-- +module Math.Grads.Angem.Internal.MatrixOperations + ( alignmentFunc + , rotation2D + ) where + +import Linear.Matrix (M22, det22, + transpose, (!*!), + (*!)) +import Linear.V2 (V2 (..)) +import Linear.Vector (negated, (^+^), + (^-^)) +import Math.Grads.Angem.Internal.VectorOperations (avg) + +-- | Given two lists of points produces function that transforms coordinates of given point +-- according to allignment of first list of points on second. +-- +alignmentFunc :: [V2 Float] -> [V2 Float] -> V2 Float -> V2 Float +alignmentFunc points1 points2 = transformFunc + where + (rotationM, transitionV) = superImpose points1 points2 + transformFunc = transform rotationM transitionV + +superImpose :: [V2 Float] -> [V2 Float] -> (M22 Float, V2 Float) +superImpose points1 points2 = (rotation, transition) + where + (avg1, moved1) = moveToCenter points1 + (avg2, moved2) = moveToCenter points2 + aMatrix = transpose moved2 !*! moved1 + (u, vt) = svd aMatrix + + rotation' = rotationMatrix vt u + rotation = if det22 rotation' >= 0 + then rotation' + else case vt of + (V2 v1 v2) -> rotationMatrix (V2 v1 (negated v2)) u + + transition = avg1 - (avg2 *! rotation) + +svd :: M22 Float -> (M22 Float, M22 Float) +svd aMatrix' = (doubleToFloatM22 rotationA, doubleToFloatM22 rotationB) + where + V2 (V2 a b) (V2 c d) = floatToDoubleM22 aMatrix' + e = (a + d) / 2 + f = (a - d) / 2 + g = (c + b) / 2 + h = (c - b) / 2 + q = sqrt (e ** 2 + h ** 2) + r = sqrt (f ** 2 + g ** 2) + a1 = atan2 g f + a2 = atan2 h e + sy = q - r + s = if sy < 0 then -1 else 1 + theta = (a2 - a1) / 2 + phi = (a2 + a1) / 2 + + rotationA = V2 (V2 (cos phi) (- s * sin phi)) (V2 (sin phi) (s * cos phi)) + rotationB = V2 (V2 (cos theta) (- sin theta)) (V2 (sin theta) (cos theta)) + +moveToCenter :: [V2 Float] -> (V2 Float, [V2 Float]) +moveToCenter points = (avgPoint, movedPoints) + where + avgPoint = avg points + movedPoints = (^-^ avgPoint) <$> points + +rotationMatrix :: M22 Float -> M22 Float -> M22 Float +rotationMatrix vt u = transpose $ transpose vt !*! transpose u + +-- | Given angle in degrees produces rotation matrix that corresponds to that angle. +-- +rotation2D :: Float -> M22 Float +rotation2D angle = V2 (V2 (cos trueAngle) (- sin trueAngle)) (V2 (sin trueAngle) (cos trueAngle)) + where + trueAngle = 2 * pi * angle / 360.0 + +transform :: M22 Float -> V2 Float -> V2 Float-> V2 Float +transform rotationM transitionV = convFunc + where + convFunc = transformVector rotationM transitionV + +transformVector :: M22 Float -> V2 Float -> V2 Float-> V2 Float +transformVector rotationM transitionV v = (v *! rotationM) ^+^ transitionV + +doubleToFloatM22 :: M22 Double -> M22 Float +doubleToFloatM22 (V2 a' b') = V2 (realToFrac <$> a') (realToFrac <$> b') + +floatToDoubleM22 :: M22 Float -> M22 Double +floatToDoubleM22 (V2 a' b') = V2 (realToFrac <$> a') (realToFrac <$> b') diff --git a/src/Math/Grads/Angem/Internal/VectorOperations.hs b/src/Math/Grads/Angem/Internal/VectorOperations.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Angem/Internal/VectorOperations.hs @@ -0,0 +1,92 @@ +-- | Some useful functions for operations with vectors. +-- +module Math.Grads.Angem.Internal.VectorOperations + ( areIntersected + , avg + , eqV2 + , reflectPoint + ) where + +import Linear.Metric (distance, norm) +import Linear.V2 (V2 (..)) +import Linear.Vector ((*^), (^+^), (^/)) + +-- | End of each line shouldn't be closer then this to other line. +-- +eps :: Float +eps = 5 + +-- | Checks whether two lines intersect. +-- +areIntersected :: (V2 Float, V2 Float) -> (V2 Float, V2 Float) -> Bool +areIntersected (x@(V2 x0 y0), y@(V2 x1 y1)) (x'@(V2 x0' y0'), y'@(V2 x1' y1')) = res + where + epsA = 20 -- Minimal distance between two lines + + a = x0 * y1 - y0 * x1 + b = x0' * y1' - x1' * y0' + + x01 = x0 - x1 + x01' = x0' - x1' + y01 = y0 - y1 + y01' = y0' - y1' + + division = x01 * y01' - y01 * x01' + + px = (a * x01' - x01 * b) / division + py = (a * y01' - y01 * b) / division + + notCommonPoint = not (x `eqV2` x' || x `eqV2` y' || y `eqV2` x' || y `eqV2` y') + + inXBounds = min x0 x1 - eps < px && px < max x0 x1 + eps && min x0' x1' - eps < px && px < max x0' x1' + eps + inYBounds = min y0 y1 - eps < py && py < max y0 y1 + eps && min y0' y1' - eps < py && py < max y0' y1' + eps + + pointOnLine = pointBelongsToLine (x', y') x || pointBelongsToLine (x', y') y + || pointBelongsToLine (x, y) x' || pointBelongsToLine (x, y) y' + notDistantEnough = not (norm (x - x') > epsA && norm (x - y') > epsA && norm (y - x') > epsA && norm (y - y') > epsA) + + res = notCommonPoint && (division /= 0 && inXBounds && inYBounds || pointOnLine || notDistantEnough) + +-- | Reflects point over given line. +-- +reflectPoint :: (V2 Float, V2 Float) -> V2 Float -> V2 Float +reflectPoint (coordA, coordB) thisPoint = res + where + V2 dirA dirB = coordB - coordA + + a' = V2 (-dirB) dirA + a = a' ^/ distance a' (V2 0.0 0.0) + b' = V2 dirB (-dirA) + b = b' ^/ distance b' (V2 0.0 0.0) + + distanceFrom = distanceFromPointToLine (coordA, coordB) + normVec = if distanceFrom (thisPoint + a) < distanceFrom (thisPoint + b) then a + else b + + transform x = x + 2 * distanceFromPointToLine (coordA, coordB) x *^ normVec + + res = if pointBelongsToLine (coordA, coordB) thisPoint then thisPoint + else transform thisPoint + +distanceFromPointToLine :: (V2 Float, V2 Float) -> V2 Float -> Float +distanceFromPointToLine (V2 x1 y1, V2 x2 y2) (V2 x0 y0) = res + where + res = abs ((y2 - y1) * x0 - (x2 - x1) * y0 + x2 * y1 - y2 * x1) / sqrt ((x1 - x2) ** 2 + (y1 - y2) ** 2) + +pointBelongsToLine :: (V2 Float, V2 Float) -> V2 Float -> Bool +pointBelongsToLine (V2 x0 y0, V2 x1 y1) (V2 x' y') = (x0 * (x' - x1) + y0 * (y' - y1)) `eqFloat` 0.0 && + (min x0 x1 < x' && x' < max x0 x1 && min y0 y1 < y' && y' < max y0 y1) + +-- | Given list of points calculates centroid of these points. +-- +avg :: [V2 Float] -> V2 Float +avg points = foldl1 (^+^) points ^/ fromIntegral (length points) + +-- | Checks two vectors of coordinates for equality. +-- +eqV2 :: V2 Float -> V2 Float -> Bool +eqV2 (V2 a b) (V2 a' b') = a `eqFloat` a' && b `eqFloat` b' + +-- TODO: We need to somehow consider length of line when comparing coordinates of two points +eqFloat :: Float -> Float -> Bool +eqFloat x y = abs (x - y) < eps diff --git a/src/Math/Grads/Drawing/Coords.hs b/src/Math/Grads/Drawing/Coords.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Coords.hs @@ -0,0 +1,74 @@ +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +-- | Module providing functions for obtaining coordinates of 'GenericGraph's on +-- a 2D plane. +-- +module Math.Grads.Drawing.Coords + ( Coord + , CoordList + , CoordMap + , Drawable (..) + , EdgeFixator + , bondLength + , getCoordsForGraph + ) where + +import Control.Monad (join) +import Data.Map.Strict (keys, + singleton) +import Math.Grads.Algo.Cycles (findCycles) +import Math.Grads.Algo.Interaction (getIndices) +import Math.Grads.Drawing.Internal.Coords (Coord, + CoordList, + CoordMap, + bondLength, + coordListForDrawing) +import Math.Grads.Drawing.Internal.Cycles (getCoordsOfGlobalCycle) +import Math.Grads.Drawing.Internal.CyclesPathsAlignment (alignCyclesAndPaths) +import Math.Grads.Drawing.Internal.Paths (findPaths, getCoordsOfPath) +import Math.Grads.Drawing.Internal.Sampling (EdgeFixator, + bestSample) +import Math.Grads.GenericGraph (GenericGraph) +import Math.Grads.Graph (EdgeList, + Graph, + toList) +import System.Random (StdGen) + +-- | Type class that defines whether graph can be drawn or not. +-- +class Graph g => Drawable g v e where + -- | Change coordinates and fixate edges that shouldn't take part in sampling. + -- + edgeFixator :: g v e -> EdgeFixator e + edgeFixator = const $ (,) [] + +-- | Given 'StdGen' returns 'CoordMap', which keys correspond to indices of +-- vertices of given 'GenericGraph'. Works only for simple planar graphs. If graph +-- is neither simple nor planar, returns Nothing. This function is best used for +-- graphs that can be represented as systems of conjugated cycles and paths between +-- them. If graph contains too complex conjugated cycles, function will return Nothing. +-- +getCoordsForGraph :: (Ord v, Ord e, Eq e, Drawable GenericGraph v e) => StdGen -> GenericGraph v e -> Maybe CoordMap +getCoordsForGraph stdGen graph = if length vertices == 1 then Just (singleton 0 (0, 0)) + else res + where + (vertices, edges) = toList graph + (globalCycles, paths) = splitIntoCyclesAndPaths edges + + globalCyclesWithCoords = sequence (fmap (getCoordsOfGlobalCycle pathsWithCoords) globalCycles) + pathsWithCoords = fmap getCoordsOfPath paths + + finalCoords = join (fmap (alignCyclesAndPaths pathsWithCoords) globalCyclesWithCoords) + resCoords = join (fmap (bestSample stdGen (edgeFixator graph) (concat paths)) finalCoords) + + resMap = fmap coordListForDrawing resCoords + + res = if fmap (length . keys) resMap == pure (length vertices) then resMap else Nothing + +splitIntoCyclesAndPaths :: (Ord e, Eq e) => EdgeList e -> ([EdgeList e], [EdgeList e]) +splitIntoCyclesAndPaths edges = (globalCycles, paths) + where + globalCycles = findCycles edges + forPaths = filter (`notElem` concat globalCycles) edges + paths = findPaths forPaths $ concatMap getIndices globalCycles diff --git a/src/Math/Grads/Drawing/Internal/Coords.hs b/src/Math/Grads/Drawing/Internal/Coords.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/Coords.hs @@ -0,0 +1,87 @@ +-- | Module providing functions for working with coordinates in Drawing module. +-- +module Math.Grads.Drawing.Internal.Coords + ( Coord + , CoordList + , CoordMap + , Link + , bondLength + , coordListForDrawing + , coordListToMap + , coordMapToCoordList + , tupleToList + ) where + +import Control.Arrow ((***)) +import Data.List (sortOn) +import Data.Map.Strict (Map, fromList, (!)) +import Linear.Metric (distance, norm) +import Linear.V2 (V2 (..)) +import Linear.Vector ((^/)) +import Math.Grads.Angem (alignmentFunc) +import Math.Grads.Drawing.Internal.Utils (Coord, CoordList, pairToV2, + tupleToList, uV2) +import Math.Grads.Graph (EdgeList, GraphEdge) + +-- | (Number of vertex, edge) for linked paths. +-- +type Link e = (Int, GraphEdge e) + +-- | Map that matches indexes of vertices to coordinates of these vertices. +-- +type CoordMap = Map Int (Float, Float) + +-- | This constant is used to determine length of one edge when graph is drawn. +-- +bondLength :: Float +bondLength = 100.0 + +-- | Given 'CoordMap' and 'EdgeList' constructs 'CoordList'. +-- +coordMapToCoordList :: CoordMap -> EdgeList e -> CoordList e +coordMapToCoordList coordMap = fmap (\bond@(a, b, _) -> (bond, (toV2Coord a, toV2Coord b))) + where + toV2Coord :: Int -> V2 Float + toV2Coord = pairToV2 . (coordMap !) + +-- | Converts 'CoordList' int 'CoordMap'. +-- +coordListForDrawing :: Eq e => CoordList e -> CoordMap +coordListForDrawing coordinates = uV2 <$> coordListToMap coordsT + where + coordsT = rotateAlongLongestDist coordinates + +-- | Converts 'CoordList' to 'Map Int (V2 Float)'. +-- +coordListToMap :: Eq e => CoordList e -> Map Int (V2 Float) +coordListToMap coordinates = fromList (helper coordinates [] []) + where + + helper :: CoordList e -> [Int] -> [(Int, V2 Float)] -> [(Int, V2 Float)] + helper [] _ res = res + helper (((a, b, _), (cA, cB)) : xs) taken res | a `elem` taken && b `elem` taken = helper xs taken res + | a `elem` taken && b `notElem` taken = helper xs (b : taken) ((b, cB) : res) + | a `notElem` taken && b `elem` taken = helper xs (a : taken) ((a, cA) : res) + | otherwise = helper xs (a : b : taken) ((a, cA) : (b, cB) : res) + +rotateAlongLongestDist :: CoordList e -> CoordList e +rotateAlongLongestDist coordinates = res + where + coordsU = getFloats coordinates + (distA, distB) = findTwoMostDistantPoints coordsU + dirVec = distB - distA + + alFunc = alignmentFunc [V2 0 0, V2 1 0] [V2 0 0, dirVec ^/ norm dirVec] + res = fmap (alFunc *** alFunc) <$> coordinates + +getFloats :: CoordList e -> [V2 Float] +getFloats coords = foldl (\x y -> x ++ tupleToList y) [] (fmap snd coords) + +findTwoMostDistantPoints :: [V2 Float] -> (V2 Float, V2 Float) +findTwoMostDistantPoints points = res + where + res = head (sortOn (\(a, b) -> -(distance a b)) (allPairs points)) + + allPairs :: [a] -> [(a, a)] + allPairs [] = [] + allPairs (x : xs) = fmap (\x' -> (x, x')) xs ++ allPairs xs diff --git a/src/Math/Grads/Drawing/Internal/Cycles.hs b/src/Math/Grads/Drawing/Internal/Cycles.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/Cycles.hs @@ -0,0 +1,222 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Module that calculates coordinates of systems of conjugated cycles in graph. +-- +module Math.Grads.Drawing.Internal.Cycles + ( getCoordsOfGlobalCycle + ) where + +import qualified Data.Array as A +import Data.List (find, + groupBy, nub, + sortOn) +import qualified Data.Map.Strict as M +import Data.Maybe (catMaybes, + fromJust, + isJust, + mapMaybe) +import Linear.Metric (distance, + norm) +import Linear.V2 (V2 (..)) +import Linear.Vector ((*^), (^/)) +import Math.Grads.Algo.Cycles (findLocalCycles) +import Math.Grads.Algo.Interaction (getEnds, + getIndices) +import Math.Grads.Algo.Paths (findBeginnings) +import Math.Grads.Algo.Traversals (dfsCycle) +import Math.Grads.Angem (alignmentFunc) +import Math.Grads.Drawing.Internal.Coords (Link, + bondLength) +import Math.Grads.Drawing.Internal.CyclesPathsAlignment (bondsToAlignTo, + bondsToAlignToExtreme) +import Math.Grads.Drawing.Internal.Utils (Coord, + CoordList, + centroid, + cleanCoordList, + cleanListOfCoordLists, + compareCoords, + findIncidentCoords, + reflectCycle, + tupleToList) +import Math.Grads.GenericGraph (gAdjacency, + gIndex) +import Math.Grads.Graph (EdgeList, + GraphEdge, + fromList, + vCount) +import Math.Grads.Utils (uniter) + +-- | Calculates coordinates of system of cycles and coordinates of edges that are adjacent to it. +-- +getCoordsOfGlobalCycle :: Eq e => [CoordList e] -> EdgeList e -> Maybe (CoordList e) +getCoordsOfGlobalCycle paths globalCycle = if not (null localCycles) && isJust alignedM then Just res + else Nothing + where + localCycles = findLocalCycles globalCycle + + localCyclesWithCoords = sortOn (\x -> - (length x)) (fmap getCoordsOfLocalCycle localCycles) + alignedM = greedyAlignmentOfLocalCycles [head localCyclesWithCoords] (tail localCyclesWithCoords) + + aligned = fromJust alignedM + cleanAligned = cleanCoordList (concat aligned) [] + + res = restoreEndsForCycle cleanAligned paths aligned + +getCoordsOfLocalCycle :: EdgeList e -> CoordList e +getCoordsOfLocalCycle thisCycle = matchBonds thisCycle (getCoordsOfPolygon (length thisCycle)) + +getCoordsOfPolygon :: Int -> [(V2 Float, V2 Float)] +getCoordsOfPolygon number = let coords = fmap getPoint [0..number - 1] in (last coords, head coords) : uniter coords + where + angle = 2 * pi / fromIntegral number + radius = bondLength / sin (angle / 2) / 2 + + getPoint :: Int -> V2 Float + getPoint step = V2 (radius * cos (fromIntegral step * angle)) (radius * sin (fromIntegral step * angle)) + +uniteLocalCyclesOnBond :: Coord e -> Coord e -> CoordList e -> CoordList e +uniteLocalCyclesOnBond (_, coords) (_, coords') toTransformCoords = transformFuncCoord <$> toTransformCoords + where + transformFunc' = alignmentFunc (tupleToList coords) (tupleToList coords') + transformFuncCoord (bond, (a, b)) = (bond, (transformFunc' a, transformFunc' b)) + +matchBonds :: EdgeList e -> [(V2 Float, V2 Float)] -> CoordList e +matchBonds bonds coords = matchBonds' bonds (zip bondsInd coords) + where + vertices = nub $ concatMap getEnds bonds + index = M.fromList (zip vertices [0..]) + graph = fromList (vertices, fmap (\(a, b, t) -> (index M.! a, index M.! b, t)) bonds) + graphArray = fmap fst <$> gAdjacency graph + + inds = (gIndex graph A.!) <$> dfsCycle graphArray [0 .. (vCount graph - 1)] [] + bondsInd = uniter inds ++ [(last inds, head inds)] + +matchBonds' :: EdgeList e -> [((Int, Int), (V2 Float, V2 Float))] -> CoordList e +matchBonds' bonds match = fmap (changeCoords match) bonds + +changeCoords :: [((Int, Int), (V2 Float, V2 Float))] -> GraphEdge e -> Coord e +changeCoords [] _ = error "No matching coords in changeCoords function." +changeCoords (((a', b'), (left, right)) : xs) bond@(a, b, _) | a == a' && b == b' = (bond, (left, right)) + | a == b' && b == a' = (bond, (right, left)) + | otherwise = changeCoords xs bond + +greedyAlignmentOfLocalCycles :: forall e. Eq e => [CoordList e] -> [CoordList e] -> Maybe [CoordList e] +greedyAlignmentOfLocalCycles mainCycles [] = Just mainCycles +greedyAlignmentOfLocalCycles mainCycles xs = if isJust idOfNeighborM then res + else Nothing + where + neighborExists = fmap checkForNeighbor xs + idOfNeighborM = helper neighborExists 0 + + idOfNeighbor = fromJust idOfNeighborM + neighbor = (xs !! idOfNeighbor) + + x = concat mainCycles + matches = catMaybes (concatMap findMatchingBond x) + (coordsA, coordsB) = head matches + + reflectedIfNeeded = reflectIfIntersects (uniteLocalCyclesOnBond coordsA coordsB neighbor) mainCycles (snd coordsA) + finalCycle = correctLeftMatches (snd <$> tail matches) reflectedIfNeeded x + + res = greedyAlignmentOfLocalCycles (finalCycle : mainCycles) (take idOfNeighbor xs ++ drop (idOfNeighbor + 1) xs) + + findMatchingBond :: Coord e -> [Maybe (Coord e, Coord e)] + findMatchingBond thisBond = fmap (hasMatch thisBond) neighbor + + hasMatch :: Coord e -> Coord e -> Maybe (Coord e, Coord e) + hasMatch thisBond otherBond = if compareCoords thisBond otherBond then Just (thisBond, otherBond) + else Nothing + + checkForNeighbor :: CoordList e -> Bool + checkForNeighbor = any (\otherCoord -> any (compareCoords otherCoord) x) + + helper :: [Bool] -> Int -> Maybe Int + helper [] _ = Nothing -- No neighbors for cycle in one conjugated cycle with it. + -- Theoretically it is impossible, but due to the nature of our findLocalCycles function this can happen + helper (y : ys) counter = if y then Just counter else helper ys (counter + 1) + +reflectIfIntersects :: CoordList e -> [CoordList e] -> (V2 Float, V2 Float) -> CoordList e +reflectIfIntersects thisCycle allCycles (coordA, coordB) = if intersects then reflectCycle thisCycle (coordA, coordB) + else thisCycle + where + thisCentroid = centroid thisCycle + otherCentroids = centroid <$> allCycles + intersects = any (\x -> distance x thisCentroid <= bondLength) otherCentroids + +correctLeftMatches :: forall e. Eq e => [Coord e] -> CoordList e -> CoordList e -> CoordList e +correctLeftMatches [] thisCycle _ = thisCycle +correctLeftMatches ((bond@(beg, end, _), _) : xs) thisCycle mainCycles = correctLeftMatches xs thisCycleUpdated mainCycles + where + thisCycleUpdated = catMaybes (fmap correctMatch thisCycle) + + correctMatch :: Coord e -> Maybe (Coord e) + correctMatch coord@(bond'@(a, b, t'), (coordA, coordB)) | bond == bond' = Nothing + | beg == a || end == a = Just ((a, b, t'), (substitute coordA a, coordB)) + | beg == b || end == b = Just ((a, b, t'), (coordA, substitute coordB b)) + | otherwise = Just coord + + substitute :: V2 Float -> Int -> V2 Float + substitute varCoord endToFix = + let + x = mapMaybe (helper endToFix) mainCycles + in if not (null x) then head x + else varCoord + + helper :: Int -> Coord e -> Maybe (V2 Float) + helper endToFix ((a', b', _), (coordA', coordB')) | a' == endToFix = Just coordA' + | b' == endToFix = Just coordB' + | otherwise = Nothing + +restoreEndsForCycle :: Eq e => CoordList e -> [CoordList e] -> [CoordList e] -> CoordList e +restoreEndsForCycle thisCycle [[]] _ = thisCycle +restoreEndsForCycle thisCycle paths localCycles = thisCycle ++ concat neighbors + where + verticesOfCycle = getIndices (fmap fst thisCycle) + cycleLinkingCoords = mapMaybe (findLinks verticesOfCycle) paths + counted = countNeighbors' cycleLinkingCoords + neighbors = fmap (getLinksWithCoords thisCycle localCycles) counted + +countNeighbors' :: [(Int, GraphEdge e)] -> [(Int, EdgeList e)] +countNeighbors' list = (\x -> let (a, b) = unzip x in (head a, b)) <$> groupBy (\a b -> fst a == fst b) list + +findLinks :: [Int] -> CoordList e -> Maybe (Link e) +findLinks verticesOfCycle path = if not (null found) then Just (foundVertex, fst (fromJust bond)) + else Nothing + where + found = filter (`elem` verticesOfCycle) (findBeginnings (fmap fst path)) + foundVertex = head found + bond = find (\((a, b, _), _) -> a == foundVertex || b == foundVertex) path + +getLinksWithCoords :: forall e. Eq e => CoordList e -> [CoordList e] -> (Int, EdgeList e) -> CoordList e +getLinksWithCoords thisCycle localCycles (ind, bonds) = res + where + found = findAdjacentBondsCycles thisCycle localCycles ind + + bondsLength = length bonds + alignedBonds = either (\(f, s) -> bondsToAlignTo f s bondsLength) (flip bondsToAlignToExtreme bondsLength) found + + res = assignCoords bonds alignedBonds ind + + assignCoords :: EdgeList e -> [(V2 Float, V2 Float)] -> Int -> CoordList e + assignCoords [] _ _ = [] + assignCoords (x@(a, _, _) : xs) (y@(left, right) : ys) start = if a == start then (x, y) : assignCoords xs ys start + else (x, (right, left)) : assignCoords xs ys start + assignCoords _ _ _ = error "Can not assign coords while restoring ends for cycle." + +findAdjacentBondsCycles :: forall e. Eq e => CoordList e -> [CoordList e] -> Int -> Either (Coord e, Coord e) (V2 Float, V2 Float) +findAdjacentBondsCycles bondsOfCycle cycles ind = if length neighbors == 2 then Left (leftNeighbor, rightNeighbor) + else Right (beginning, beginning + bondLength *^ direction ^/ norm direction) + where + neighbors = findIncidentCoords ind bondsOfCycle + + [leftNeighbor, rightNeighbor] = take 2 neighbors + + cyclesInPlay = cleanListOfCoordLists (filter (\x -> any (`elem` x) neighbors) cycles) [] + beginning = findCommonVertexCoord leftNeighbor rightNeighbor + direction = (beginning - centroid (head cyclesInPlay)) + (beginning - centroid (last cyclesInPlay)) + + findCommonVertexCoord :: Coord e -> Coord e -> V2 Float + findCommonVertexCoord ((a, b, _), (coordA, coordB)) ((a', b', _), _) | a == a' = coordA + | a == b' = coordA + | b == a' = coordB + | otherwise = coordB diff --git a/src/Math/Grads/Drawing/Internal/CyclesPathsAlignment.hs b/src/Math/Grads/Drawing/Internal/CyclesPathsAlignment.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/CyclesPathsAlignment.hs @@ -0,0 +1,137 @@ +-- | Module that is responsible for linking systems of conjugated cycles in graph +-- with paths between them. +-- +module Math.Grads.Drawing.Internal.CyclesPathsAlignment + ( alignCyclesAndPaths + , bondsToAlignTo + , bondsToAlignToExtreme + ) where + +import Control.Arrow (first, second, (***)) +import Data.Either (partitionEithers) +import Data.List (find) +import Data.Maybe (catMaybes, listToMaybe) +import Linear.Matrix ((!*)) +import Linear.Metric (dot, norm) +import Linear.V2 (V2 (..)) +import Linear.Vector (negated, (*^)) +import Math.Grads.Algo.Paths (findBeginnings) +import Math.Grads.Angem (alignmentFunc, rotation2D) +import Math.Grads.Drawing.Internal.Coords (bondLength) +import Math.Grads.Drawing.Internal.Utils (Coord, CoordList, + cleanCoordList, + tupleToList) +import Math.Grads.Graph (EdgeList) + +type CoordsEnds e = (CoordList e, EdgeList e) + +-- | Given cycles and paths between them unites everything into one structure if possible. +-- +alignCyclesAndPaths :: Eq e => [CoordList e] -> [CoordList e] -> Maybe (CoordList e) +alignCyclesAndPaths paths cycles = greedyAlignmentOfCyclesAndPaths (cyclesWithRestoredEnds ++ pathsWithRestoredEnds) + where + cyclesWithRestoredEnds = fmap linksForCycle cycles + pathsWithRestoredEnds = fmap linksForPath paths + + linksForCycle :: CoordList e -> (CoordList e, EdgeList e) + linksForCycle thisCycle = (thisCycle, findBondsToFind (fmap fst thisCycle)) + + linksForPath :: CoordList e -> (CoordList e, EdgeList e) + linksForPath thisPath = (thisPath, helper' (fmap fst thisPath)) + + helper' :: EdgeList e -> EdgeList e + helper' pathBondList = if length pathBondList == 1 then pathBondList + else findBondsToFind pathBondList + +greedyAlignmentOfCyclesAndPaths :: Eq e => [(CoordList e, EdgeList e)] -> Maybe (CoordList e) +greedyAlignmentOfCyclesAndPaths [] = Nothing +greedyAlignmentOfCyclesAndPaths [x] = Just (fst x) +greedyAlignmentOfCyclesAndPaths (thisPart : otherParts) = if not (null toAdd) then res + else Nothing + where + theseCoords = fst thisPart + bondsToFind = snd thisPart + alignedNeighbors = fmap (detectAndAlignNeighbors bondsToFind theseCoords) otherParts + + (toAdd, leftParts) = first concat (partitionEithers alignedNeighbors) + + newTheseCoords = cleanCoordList (toAdd ++ theseCoords) [] + + edgeList = fmap fst newTheseCoords + newBondsToFind = findBondsToFind edgeList + + res = greedyAlignmentOfCyclesAndPaths ((newTheseCoords, newBondsToFind) : leftParts) + +detectAndAlignNeighbors :: Eq e => EdgeList e -> CoordList e -> CoordsEnds e -> Either (CoordList e) (CoordsEnds e) +detectAndAlignNeighbors bondsToFind theseCoords theseCoordsEnds = maybe (Right theseCoordsEnds) Left neighsOrLeft + where + neighsOrLeft = detectAndAlignNeighborsM bondsToFind theseCoords theseCoordsEnds + +detectAndAlignNeighborsM :: Eq e => EdgeList e -> CoordList e -> CoordsEnds e -> Maybe (CoordList e) +detectAndAlignNeighborsM bondsToFind theseCoords (coords, ends) = do + let found' = catMaybes (fmap (\x -> find (== x) bondsToFind) ends) + found <- listToMaybe found' + + let findBondToAlign = find (\(a, _) -> a == found) + + alignCoords <- coordToList <$> findBondToAlign theseCoords + toAlignCoords <- coordToList <$> findBondToAlign coords + + let alignFunc = alignmentFunc alignCoords toAlignCoords + + Just (fmap (second (alignFunc *** alignFunc)) coords) + where + coordToList :: Coord e -> [V2 Float] + coordToList = tupleToList . snd + +findBondsToFind :: EdgeList e -> EdgeList e +findBondsToFind bonds = catMaybes ((\ind -> find (\(a, b, _) -> a == ind || b == ind) bonds) <$> findBeginnings bonds) + +-- | Constructs edge that will be used to align to cycle containing given 'Coord's. +-- +bondsToAlignTo :: Coord e -> Coord e -> Int -> [(V2 Float, V2 Float)] +bondsToAlignTo ((a, b, _), (pointA, pointB)) ((a', b', _), (pointA', pointB')) number = resultingVectors + where + coordA = pointB - pointA + coordB = pointB' - pointA' + ((vecA, vecB), linkingPoint) | a == a' = ((negated coordA, negated coordB), pointA) + | a == b' = ((negated coordA, coordB), pointA) + | b == a' = ((coordA, negated coordB), pointB) + | otherwise = ((coordA, coordB), pointB) + + direction' = vecA + vecB + direction = (bondLength / norm direction') *^ direction' + toTopAngle = (180.0 - 180.0 * acos (dot vecA vecB / (norm vecA * norm vecB)) / pi) / 2.0 + angle' = 180.0 / fromIntegral number + startingAngle = (180.0 - (fromIntegral number - 1.0) * angle') / 2.0 + + dirA = dot (start (toTopAngle + startingAngle)) direction + dirB = dot (start (-(toTopAngle + startingAngle))) direction + startingPoint | dirA >= 0 && dirB >= 0 && dirA > dirB = start (toTopAngle + startingAngle) + | dirA >= 0 && dirB >= 0 = start (-(toTopAngle + startingAngle)) + | dirA >= 0 = start (toTopAngle + startingAngle) + | otherwise = start (-(toTopAngle + startingAngle)) + + mult = if dot (start (toTopAngle + startingAngle)) direction > 0 then 1 else (-1) + resultingVectors = (\x -> (linkingPoint, linkingPoint + x)) <$> getDirections startingPoint 1 angle' number mult + + start :: Float -> V2 Float + start angle = rotation2D angle !* ((bondLength / norm vecA) *^ negated vecA) + +-- | If we have complicated situation where we need to calculate bonds to align to +-- for vertex in cycle that has more then 2 neighbors then we pass direction in +-- which we want to place neighbors and use bondsToAlignToExtreme function. +-- Otherwise we use bondsToAlignTo function. +-- +bondsToAlignToExtreme :: (V2 Float, V2 Float) -> Int -> [(V2 Float, V2 Float)] +bondsToAlignToExtreme (beg, end) number = resultingVectors + where + direction = end - beg + startingPointComplicated = rotation2D (-40.0) !* ((bondLength / norm direction) *^ direction) + resultingVectors = (\x -> (beg, beg + x)) <$> getDirections startingPointComplicated 1 47.0 number 1 + +getDirections :: V2 Float -> Int -> Float -> Int -> Float -> [V2 Float] +getDirections prev counter angle number mult = if counter < number then prev : getDirections new (counter + 1) angle number mult + else [prev] + where + new = rotation2D (mult * angle) !* prev diff --git a/src/Math/Grads/Drawing/Internal/Paths.hs b/src/Math/Grads/Drawing/Internal/Paths.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/Paths.hs @@ -0,0 +1,286 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Module that calculates coordinates of paths between systems of conjugated cycles in graph. +-- +module Math.Grads.Drawing.Internal.Paths + ( findPaths + , getCoordsOfPath + ) where + +import Data.List (delete, find, + intersect, + maximumBy, + union, (\\)) +import Data.Map.Strict (fromList, + (!)) +import Data.Maybe (catMaybes, + fromJust, + isJust) +import Data.Ord (comparing) +import Linear.V2 (V2 (..)) +import Math.Grads.Algo.Interaction (getIndices, getSharedVertex, + getVertexIncident, + isIncident) +import Math.Grads.Algo.Paths (dfsSearch, findBeginnings) +import Math.Grads.Algo.Traversals (dfs) +import Math.Grads.Angem (alignmentFunc) +import Math.Grads.Drawing.Internal.Coords (Link, + bondLength) +import Math.Grads.Drawing.Internal.CyclesPathsAlignment (bondsToAlignTo) +import Math.Grads.Drawing.Internal.Utils (Coord, + CoordList, + findIncidentCoords, + tupleToList) +import Math.Grads.Graph (EdgeList, + GraphEdge) +type BondV2 = (V2 Float, V2 Float) + +type PathWithLinks e = (CoordList e, [Link e]) + +data Path e = Path { + pStart :: Int, + pFinish :: Int, + pBonds :: EdgeList e +} + +-- | Calculates coordinates of path. +-- +getCoordsOfPath :: forall e. Eq e => EdgeList e -> CoordList e +getCoordsOfPath bonds = fst (greedy pathWithCoords) + where + paths = splitPathIntoLongest bonds [] + pathsWithCoords = fmap (\Path { pBonds=bonds' } -> pathToCoords bonds') paths + pathWithCoords = zip pathsWithCoords (getLinks paths []) + + getLinks :: [Path e] -> [Int] -> [[Link e]] + getLinks (x : xs) taken = + let + links = findLinkingPoints x taken xs + in links : getLinks xs (taken `union` fmap fst (countNeighbors links)) + getLinks _ _ = error "No links for path." + + greedy :: [PathWithLinks e] -> PathWithLinks e + greedy [x] = x + greedy (x : xs) = greedy (uniteOnLinkingBonds x xs : filter (helper' x) xs) + greedy _ = error "Greedy function on an empty list." + + helper' :: PathWithLinks e -> PathWithLinks e -> Bool + helper' x' y' = null $ (snd <$> snd x') `intersect` (fst <$> fst y') + +-- | Takes path, list of vertices which have been processed and returns links for each path. +-- +findLinkingPoints :: forall e. Path e -> [Int] -> [Path e] -> [Link e] +findLinkingPoints Path { pBonds=list } taken = helper + where + getInc :: Int -> EdgeList e + getInc n = filter (`isIncident` n) list + + couldTake :: Int -> Bool + couldTake n = n `notElem` taken && (not . null $ getInc n) + + helper :: [Path e] -> [Link e] + helper [] = [] + helper (Path beg end list' : xs) | couldTake beg = wrapResult beg list' : helper xs + | couldTake end = wrapResult end list' : helper xs + | otherwise = helper xs + + wrapResult :: Int -> EdgeList e -> Link e + wrapResult n l = (fromJust $ getSharedVertex b1 b2, b2) + where + b1 = head $ getInc n + b2 = head $ getVertexIncident l n + +splitPathIntoLongest :: Eq e => EdgeList e -> [Int] -> [Path e] +splitPathIntoLongest [] _ = [] +splitPathIntoLongest bonds taken = firstPath : splitPathIntoLongest restBonds newTaken + where + ends' = findBeginnings bonds + ends = filter (`notElem` taken) ends' + + startEnds = if not (null taken) then taken else ends' + + allPaths = filter (not . null) ((\(x, y) -> maybe [] fst (dfsSearch bonds x y)) <$> allPairs startEnds ends) + allPathsTrue = concatMap (\x -> maybe [x] (splitOnPoint x) (findPointsToSplit x taken)) allPaths + + longestPath = maximumBy (comparing length) allPathsTrue + [start, finish] = findBeginnings longestPath + + restBonds = bonds \\ longestPath + firstPath = Path { pStart=start, pFinish=finish, pBonds=longestPath } + newTaken = taken `union` getIndices longestPath + +pathToCoords :: forall e. EdgeList e -> CoordList e +pathToCoords bonds = matchBondsOfPath coords bonds + where + angle = pi / 6.0 + radius = bondLength + allPoints = getPoint <$> concat (repeat [1.0, -1.0]) + dfsRes = dfs bonds (head (findBeginnings bonds)) + coords = zip (getIndicesEdges dfsRes) (buildPath allPoints) + + getPoint :: Float -> V2 Float + getPoint m = V2 (radius * cos (m * angle)) (radius * sin (m * angle)) + + getIndicesEdges :: EdgeList e -> [Int] + getIndicesEdges [] = error "Get indices edges on empy list." + getIndicesEdges [(a, b, _)] = [a, b] + getIndicesEdges (bnd@(a, b, _) : bnd'@(a', b', _) : xs) = if a == a' || a == b' then b : a : helper (bnd' : xs) bnd + else a : b : helper (bnd' : xs) bnd + + helper :: EdgeList e -> GraphEdge e -> [Int] + helper [] _ = error "Get indices edges helper on empty list." + helper [(a', b', _)] (a, b, _) = if a' == a || a' == b then [b'] + else [a'] + helper (bnd@(a, b, _) : bnd'@(a', b', _) : xs) _ = if a == a' || a == b' then a : helper (bnd' : xs) bnd + else b : helper (bnd' : xs) bnd + +buildPath :: [V2 Float] -> [V2 Float] +buildPath [] = [] +buildPath (y : ys) = V2 0.0 0.0 : y : helper ys y + where + helper :: [V2 Float] -> V2 Float -> [V2 Float] + helper [] _ = [] + helper (b:xs) b' = let newCoords = b + b' in newCoords : helper xs newCoords + +matchBondsOfPath :: forall e. [(Int, V2 Float)] -> EdgeList e -> CoordList e +matchBondsOfPath matches = helper + where + mapOfCoords = fromList matches + + helper :: EdgeList e -> CoordList e + helper [] = [] + helper (bond@(a, b, _) : xs) = (bond, (mapOfCoords ! a, mapOfCoords ! b)) : helper xs + +uniteOnLinkingBonds :: forall e. PathWithLinks e -> [PathWithLinks e] -> PathWithLinks e +uniteOnLinkingBonds (mainPath, uniteThis) otherPaths = pathUniter coordsToAdd (mainPath, []) + where + counted = countNeighbors uniteThis + neighbors = fmap getCoordsOfLinks counted + + coordsToAdd = concatMap align neighbors + + align :: (Int, [BondV2]) -> [PathWithLinks e] + align (ind, toAlignBonds) = alignOnBonds <$> zip toAlignBonds (findNeighbors otherPaths ind) + + getCoordsOfLinks :: (Int, Int) -> (Int, [BondV2]) + getCoordsOfLinks (ind, counts) = + let + (leftBond, rightBond) = findAdjacent mainPath ind + in (ind, bondsToAlignTo leftBond rightBond counts) + + pathUniter :: [PathWithLinks e] -> PathWithLinks e -> PathWithLinks e + pathUniter [] res = res + pathUniter ((a, b) : xs) (resA, resB) = pathUniter xs (resA ++ a, resB ++ b) + +countNeighbors :: forall e. [Link e] -> [(Int, Int)] +countNeighbors list = zip allInds (fmap (numberOfNeighbors list 0) allInds) + where + allInds = linkingIndices list [] + + linkingIndices :: [Link e] -> [Int] -> [Int] + linkingIndices [] res = res + linkingIndices (x : xs) res = if fst x `elem` res then linkingIndices xs res + else linkingIndices xs (fst x : res) + + numberOfNeighbors :: [Link e] -> Int -> Int -> Int + numberOfNeighbors [] acc _ = acc + numberOfNeighbors (x : xs) acc ind = if fst x == ind then numberOfNeighbors xs (acc + 1) ind + else numberOfNeighbors xs acc ind + +splitOnMultiplePoints :: forall e. Eq e => EdgeList e -> [Int] -> [EdgeList e] +splitOnMultiplePoints bonds cut = helper [bonds] cut [] + where + helper :: [EdgeList e] -> [Int] -> [EdgeList e] -> [EdgeList e] + helper [] _ _ = error "Split on multiple points helper on empty list." + helper lastCut [] res = res ++ lastCut + helper (x : xs) (y : ys) res = if y `elem` getIndices x then helper (splitOnPoint x y ++ xs) ys res + else helper xs (y : ys) (x : res) + +splitOnPoint :: forall e. Eq e => EdgeList e -> Int -> [EdgeList e] +splitOnPoint list point = filter (not . null) foundNeighbors + where + foundNeighbors = fmap splitter list + + splitter :: GraphEdge e -> EdgeList e + splitter bond@(a, b, _) | a == point = bond : dfs (delete bond list) b + | b == point = bond : dfs (delete bond list) a + | otherwise = [] + +allPairs :: [Int] -> [Int] -> [(Int, Int)] +allPairs starts ends = concatMap (\start -> fmap (\end -> (start, end)) ends) starts + +-- This function is used for splitting one path into substantial pieces during calculation of coordinates of this path +findPointsToSplit :: forall e. EdgeList e -> [Int] -> Maybe Int +findPointsToSplit _ [] = Nothing +findPointsToSplit [] _ = Nothing +findPointsToSplit [_] _ = Nothing +findPointsToSplit list taken = helper ((tail . init) list) + where + helper :: EdgeList e -> Maybe Int + helper [] = Nothing + helper ((a, b, _) : xs) | a `elem` taken = Just a + | b `elem` taken = Just b + | otherwise = helper xs + +-- This function is used for splitting one path into several paths if one vertex of path belongs to cycle +findPointsToSplitHard :: forall e. Eq e => EdgeList e -> [Int] -> [Int] +findPointsToSplitHard _ [] = [] +findPointsToSplitHard [] _ = [] +findPointsToSplitHard [_] _ = [] +findPointsToSplitHard list taken = helper list [] + where + helper :: EdgeList e -> [Int] -> [Int] + helper [] acc = acc + helper (x@(a, b, _) : xs) acc | a `notElem` acc && a `elem` taken && hasNeighbor a x = helper xs (a : acc) + | b `notElem` acc && b `elem` taken && hasNeighbor b x = helper xs (b : acc) + | otherwise = helper xs acc + + hasNeighbor :: Int -> GraphEdge e -> Bool + hasNeighbor ind x = isJust (find (\bond -> bond /= x && isIncident bond ind) list) + +-- | Finds all paths between cycles in graph. +-- +findPaths :: Eq e => EdgeList e -> [Int] -> [EdgeList e] +findPaths [] _ = [] +findPaths bonds [] = [bonds] +findPaths bonds taken = allPathsTrue + where + paths = findPaths' bonds + allPathsTrue = concatMap (\x -> let cut = findPointsToSplitHard x taken in splitOnMultiplePoints x cut) paths + +findPaths' :: Eq e => EdgeList e -> [EdgeList e] +findPaths' [] = [] +findPaths' bonds@(x : _) = newPath : findPaths' (filter (`notElem` newPath) bonds) + where + newPath = findPath bonds bonds [x] x + +findPath :: Eq e => EdgeList e -> EdgeList e -> EdgeList e -> GraphEdge e -> EdgeList e +findPath [] _ found _ = found +findPath (x@(a, b, _) : xs) bonds found pathFrom@(src, dst, _) = if cond then findPath xs bonds newFound pathFrom + else findPath xs bonds found pathFrom + where + cond = (a == dst || b == src || a == src || b == dst) && (x `notElem` found) + newFound = findPath bonds bonds (found ++ [x]) x + +alignOnBonds :: (BondV2, (BondV2, PathWithLinks e)) -> PathWithLinks e +alignOnBonds (coordsA, (coordsB, (list, toSave))) = (resCoords, toSave) + where + align = alignmentFunc (tupleToList coordsA) (tupleToList coordsB) + resCoords = fmap (\(bond, (coordA', coordB')) -> (bond, (align coordA', align coordB'))) list + +findNeighbors :: [PathWithLinks e] -> Int -> [(BondV2, PathWithLinks e)] +findNeighbors [] _ = [] +findNeighbors (x : xs) ind = if not (null found) then (head found, x) : findNeighbors xs ind + else findNeighbors xs ind + where + found = catMaybes (fmap coordsOfAdjacentBond (fst x)) + + coordsOfAdjacentBond :: Coord e -> Maybe BondV2 + coordsOfAdjacentBond ((a, b, _), (coordsA, coordsB)) | a == ind = Just (coordsA, coordsB) + | b == ind = Just (coordsB, coordsA) + | otherwise = Nothing + +findAdjacent :: CoordList e -> Int -> (Coord e, Coord e) +findAdjacent list ind = (leftNeighbor, rightNeighbor) + where + [leftNeighbor, rightNeighbor] = take 2 (findIncidentCoords ind list) diff --git a/src/Math/Grads/Drawing/Internal/Sampling.hs b/src/Math/Grads/Drawing/Internal/Sampling.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/Sampling.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} + +-- | Module that provides function for sampling of coords of graph. +-- +module Math.Grads.Drawing.Internal.Sampling + ( EdgeFixator + , bestSample + ) where + +import Data.List (delete, find, (\\)) +import Data.Maybe (fromJust) +import Math.Grads.Algo.Traversals (dfs) +import Math.Grads.Angem (areIntersected, eqV2) +import Math.Grads.Drawing.Internal.Coords (CoordMap, + coordListForDrawing, + coordMapToCoordList) +import Math.Grads.Drawing.Internal.Utils (Coord, CoordList, + randomVectors, reflectBond) +import Math.Grads.Graph (EdgeList, GraphEdge) +import System.Random (StdGen) + +-- | Type alias for function that, given 'CoordMap' of graph, returns modified +-- version of that 'CoordMap' alongside with 'EdgeList' of graph's edges of that +-- shouldn't participate in sampling. +-- +type EdgeFixator e = CoordMap -> (EdgeList e, CoordMap) + +-- | Finds conformation with minimal number of intersections. +-- +bestSample :: Eq e => StdGen -> EdgeFixator e -> EdgeList e -> CoordList e -> Maybe (CoordList e) +bestSample stdGen edgeFixator bondsOfPaths coords = res + where + (fixedBonds, coordsChangedMap) = edgeFixator (coordListForDrawing coords) + + coordsChanged = coordMapToCoordList coordsChangedMap (fmap fst coords) + samples = generateSamples stdGen coordsChanged (bondsOfPaths \\ fixedBonds) + curInt = findIntersections (head samples) + + resSample = if curInt == 0 then head samples + else minInterSample (tail samples) (head samples) curInt + + res = if findIntersections resSample /= 0 then Nothing + else Just resSample + +minInterSample :: Eq e => [CoordList e] -> CoordList e -> Int -> CoordList e +minInterSample [] prev _ = prev +minInterSample (x : xs) prev prevMin | curInt' >= prevMin = minInterSample xs prev prevMin + | curInt' == 0 = x + | otherwise = minInterSample xs x curInt' + where + curInt' = findIntersections x + +generateSamples :: Eq e => StdGen -> CoordList e -> EdgeList e -> [CoordList e] +generateSamples _ coords [] = [coords] +generateSamples stdGen coords rotatableBonds = (rotateOnBonds coords <$>) filteredSubsets + where + numberOfSamples = 2000 + lengthOfBonds = length rotatableBonds + + vectors = replicate lengthOfBonds 0 : randomVectors stdGen lengthOfBonds numberOfSamples + filteredSubsets = fmap (\vector -> concatMap (\(x, y) -> [y | x == 1]) (zip vector rotatableBonds)) vectors + +rotateOnBonds :: Eq e => CoordList e -> EdgeList e -> CoordList e +rotateOnBonds = foldl rotateOnBond + +rotateOnBond :: Eq e => CoordList e -> GraphEdge e -> CoordList e +rotateOnBond coords bond = res + where + bondItself@((_, b, _), (coordA, coordB)) = fromJust (find ((== bond) . fst) coords) + + toTheRightBonds = dfs (fst <$> delete bondItself coords) b + + toTheRightCoords = filter (\(x, _) -> x `elem` toTheRightBonds) coords + toTheLeftCoords = filter (\(x, _) -> notElem x toTheRightBonds) coords + + (doNotRotate, rotate) = if length toTheLeftCoords < length toTheRightCoords then (toTheRightCoords, toTheLeftCoords) + else (toTheLeftCoords, toTheRightCoords) + + res = if null toTheLeftCoords || null toTheRightCoords then coords + else doNotRotate ++ ((`reflectBond` (coordA, coordB)) <$> rotate) + +doOverlap :: Coord e -> Coord e -> Bool +doOverlap ((a, b, _), (coordA, coordB)) ((a', b', _), (coordA', coordB')) = condA || condB + where + condA = coordA `eqV2` coordA' && coordB `eqV2` coordB' || + coordA `eqV2` coordB' && coordB `eqV2` coordA' + condB = a /= a' && coordA `eqV2` coordA' || a /= b' && coordA `eqV2` coordB' || + b /= a' && coordB `eqV2` coordA' || b /= b' && coordB `eqV2` coordB' + +findIntersections :: forall e. Eq e => CoordList e -> Int +findIntersections [] = error "Find intersections helper on empty list." +findIntersections [_] = 0 +findIntersections (x : xs) = foldl (allLeftIntersections x) 0 xs + findIntersections xs + +allLeftIntersections :: Eq e => Coord e -> Int -> Coord e -> Int +allLeftIntersections coord x' coord' = x' + addIfIntersect coord coord' + +addIfIntersect :: Eq e => Coord e -> Coord e -> Int +addIfIntersect x@(bond, coords) coord@(bond', coords') = fromEnum cond + where + cond = bond /= bond' && (doOverlap x coord || areIntersected coords coords') diff --git a/src/Math/Grads/Drawing/Internal/Utils.hs b/src/Math/Grads/Drawing/Internal/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Drawing/Internal/Utils.hs @@ -0,0 +1,121 @@ +-- | Module that provides utility functions for whole Drawing module. +-- +module Math.Grads.Drawing.Internal.Utils + ( Coord + , CoordList + , randomVectors + , findIncidentCoords + , reflectCycle + , reflectBond + , centroid + , tupleToList + , compareCoords + , cleanListOfCoordLists + , cleanCoordList + , coordElem + , pairToV2 + , uV2 + ) where + +import Control.Arrow ((***)) +import Data.List (unfoldr) +import Linear.V2 (V2 (..)) +import Math.Grads.Algo.Interaction (isIncident) +import Math.Grads.Angem (reflectPoint) +import Math.Grads.Graph (GraphEdge) +import System.Random (StdGen, randomR) + +-- | Alias for list of 'Coord's. +-- +type CoordList e = [Coord e] + +-- | 'Coord' is pair that containts graph's edge and coordinates of vertices that +-- are incident to it. +-- +type Coord e = (GraphEdge e, (V2 Float, V2 Float)) + +-- | Generates vector of random 'Int's of given length. +-- +randomVector :: Int -> StdGen -> ([Int], StdGen) +randomVector len gen = helper 0 ([], gen) + where + helper :: Int -> ([Int], StdGen) -> ([Int], StdGen) + helper currentLength (a, g) | currentLength < len = helper (currentLength + 1) (headA : a, newGen) + | otherwise = (a, g) + where + (headA, newGen) = randomR (0, 1) g + +-- | Generates list of random vectors of length numberOfVectors. +-- Each vector has length lengthOfVector. +-- +randomVectors :: StdGen -> Int -> Int -> [[Int]] +randomVectors gen lengthOfVector numberOfVectors = helper gen + where + helper = take numberOfVectors . unfoldr (Just . randomVector lengthOfVector) + +-- | Find all 'Coord's in 'CoordList' that are incident to vertex with given index. +-- +findIncidentCoords :: Int -> CoordList e -> CoordList e +findIncidentCoords ind = filter (flip isIncident ind . fst) + +-- | Reflect given cycle in form of 'CoordList' over given axis. +-- +reflectCycle :: CoordList e -> (V2 Float, V2 Float) -> CoordList e +reflectCycle thisCycle = (<$> thisCycle) . flip reflectBond + +-- | Reflect given 'Coord' over given axis. +-- +reflectBond :: Coord e -> (V2 Float, V2 Float) -> Coord e +reflectBond coord ends = fmap (reflectPoint ends *** reflectPoint ends) coord + +-- | Calculates centroid of vertices in given 'CoordList'. +-- +centroid :: CoordList e -> V2 Float +centroid coords' = sum coords / fromIntegral (length coords) + where + coords = snd <$> foldl (\x ((a, b, _), (coordA, coordB)) -> helper x (a, coordA) (b, coordB)) [] coords' + helper list (a, coordA') (b, coordB') | a `elem` fmap fst list && b `elem` fmap fst list = list + | a `elem` fmap fst list = (b, coordB') : list + | b `elem` fmap fst list = (a, coordA') : list + | otherwise = (a, coordA') : ((b, coordB') : list) + +-- | Converts tuple to lis. +-- +tupleToList :: (a, a) -> [a] +tupleToList (x, y) = [x, y] + +-- | Converts 'V2' to pair. +-- +uV2 :: V2 Float -> (Float, Float) +uV2 (V2 a b) = (a, b) + +-- | Converts pair to 'V2'. +-- +pairToV2 :: (Float, Float) -> V2 Float +pairToV2 (a, b) = V2 a b + +-- | Given 'CoordList' of conjugated cycles leaves only cycles that don't intersect +-- with each other excluding first cycle in list that is taken by default. +-- +cleanListOfCoordLists :: Eq e => [CoordList e] -> [CoordList e] -> [CoordList e] +cleanListOfCoordLists [] final = final +cleanListOfCoordLists (x : xs) [] = cleanListOfCoordLists xs [x] +cleanListOfCoordLists (x : xs) final = if any (\thisCycle -> any (`coordElem` thisCycle) x) xs then cleanListOfCoordLists xs final + else cleanListOfCoordLists xs (x : final) + +-- | Leaves only unique 'Coord's in given 'CoordList'. +-- +cleanCoordList :: Eq e => CoordList e -> CoordList e -> CoordList e +cleanCoordList [] coords = coords +cleanCoordList (x : xs) coords = if not (x `coordElem` coords) then cleanCoordList xs (x : coords) + else cleanCoordList xs coords + +-- | Checks that 'Coord' is present in 'CoordList'. +-- +coordElem :: Eq e => Coord e -> CoordList e -> Bool +coordElem coord = any (compareCoords coord) + +-- | Comparator for 'Coord's. +-- +compareCoords :: Eq e => Coord e -> Coord e -> Bool +compareCoords (a, _) (b, _) = a == b diff --git a/src/Math/Grads/GenericGraph.hs b/src/Math/Grads/GenericGraph.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/GenericGraph.hs @@ -0,0 +1,243 @@ +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE InstanceSigs #-} +{-# LANGUAGE ViewPatterns #-} + +-- | Module that provides abstract implementation of graph-like data structure +-- 'GenericGraph' and many helpful functions for interaction with 'GenericGraph'. +-- +module Math.Grads.GenericGraph + ( GenericGraph (..) + , addEdges + , addVertices + , applyG + , applyV + , getVertices + , getEdge + , isConnected + , removeEdges + , removeVertices + , safeAt + , safeIdx + , subgraph + , sumGraphs + , typeOfEdge + ) where + +import Control.Arrow (first) +import Data.Aeson (FromJSON (..), ToJSON (..), defaultOptions, + genericParseJSON, genericToJSON) +import Data.Array (Array) +import qualified Data.Array as A +import Data.List (find, groupBy, sortBy) +import Data.Map.Strict (Map, mapKeys, member, (!)) +import qualified Data.Map.Strict as M +import Data.Maybe (fromJust, fromMaybe, isJust) +import qualified Data.Set as S +import GHC.Generics (Generic) +import Math.Grads.Graph (Graph (..)) + +-- | Generic undirected graph which stores elements of type v in its vertices (e.g. labels, atoms, states etc) +-- and elements of type e in its edges (e.g. weights, bond types, functions over states etc). +-- Note that loops and multiple edges between two vertices are allowed. +-- +data GenericGraph v e = GenericGraph { gIndex :: Array Int v -- ^ 'Array' that contains vrtices of graph + , gRevIndex :: Map v Int -- ^ 'Map' that maps vertices to their indices + , gAdjacency :: Array Int [(Int, e)] -- ^ adjacency 'Array' of graph + } + deriving (Generic) + +instance (Ord v, Eq e, ToJSON v, ToJSON e) => ToJSON (GenericGraph v e) where + toJSON (toList -> l) = genericToJSON defaultOptions l + +instance (Ord v, Eq e, FromJSON v, FromJSON e) => FromJSON (GenericGraph v e) where + parseJSON v = fromList <$> genericParseJSON defaultOptions v + +instance Graph GenericGraph where + fromList :: (Ord v, Eq v) => ([v], [(Int, Int, e)]) -> GenericGraph v e + fromList (vertices, edges) = GenericGraph idxArr revMap adjArr + where + count = length vertices + idxArr = A.listArray (0, count - 1) vertices + revMap = M.fromList $ zip vertices [0..] + indices = concatMap insertFunc edges + insertFunc (at, other, b) | at == other = [(at, (other, b))] + | otherwise = [(at, (other, b)), (other, (at, b))] + adjArr = A.accumArray (flip (:)) [] (0, count - 1) indices + + toList :: (Ord v, Eq v) => GenericGraph v e -> ([v], [(Int, Int, e)]) + toList (GenericGraph idxArr _ adjArr) = (snd <$> A.assocs idxArr, edges) + where + edges = distinct . concatMap toEdges . A.assocs $ adjArr + toEdges (k, v) = map (toAscending k) v + toAscending k (a, b) | k > a = (a, k, b) + | otherwise = (k, a, b) + compare3 (at1, other1, _) (at2, other2, _) = compare (at1, other1) (at2, other2) + eq3 v1 v2 = compare3 v1 v2 == EQ + distinct = map head . groupBy eq3 . sortBy compare3 + + vCount :: GenericGraph v e -> Int + vCount (GenericGraph idxArr _ _) = length idxArr + + (!>) :: (Ord v, Eq v) => GenericGraph v e -> v -> [(v, e)] + (GenericGraph idxArr revMap adjArr) !> at = first (idxArr A.!) <$> adjacent + where + idx = revMap ! at + adjacent = adjArr A.! idx + + (?>) :: (Ord v, Eq v) => GenericGraph v e -> v -> Maybe [(v, e)] + gr@(GenericGraph _ revMap _) ?> at | at `member` revMap = Just (gr !> at) + | otherwise = Nothing + + + (!.) :: GenericGraph v e -> Int -> [(Int, e)] + (!.) (GenericGraph _ _ adjArr) = (adjArr A.!) + + (?.) :: GenericGraph v e -> Int -> Maybe [(Int, e)] + gr@(GenericGraph _ _ adjArr) ?. idx | idx `inBounds` A.bounds adjArr = Just (gr !. idx) + | otherwise = Nothing + where + -- | Check whether or not given value is betwen bounds. + -- + inBounds :: Ord a => a -> (a, a) -> Bool + inBounds i (lo, hi) = (i >= lo) && (i <= hi) + + +instance (Ord v, Eq v, Show v, Show e) => Show (GenericGraph v e) where + show gr = unlines . map fancyShow . filter (\(a, b, _) -> a < b) . snd . toList $ gr + where + idxArr = gIndex gr + fancyShow (at, other, bond) = concat [show $ idxArr A.! at, "\t", show bond, "\t", show $ idxArr A.! other] + +instance Functor (GenericGraph v) where + fmap f (GenericGraph idxArr revMap adjArr) = GenericGraph idxArr revMap (((f <$>) <$>) <$> adjArr) + + +-- | 'fmap' which acts on adjacency lists of each vertex. +-- +applyG :: ([(Int, e1)] -> [(Int, e2)]) -> GenericGraph v e1 -> GenericGraph v e2 +applyG f (GenericGraph idxArr revMap adjArr) = GenericGraph idxArr revMap (f <$> adjArr) + +-- | 'fmap' which acts on vertices. +-- +applyV :: Ord v2 => (v1 -> v2) -> GenericGraph v1 e -> GenericGraph v2 e +applyV f (GenericGraph idxArr revMap adjArr) = GenericGraph (f <$> idxArr) (mapKeys f revMap) adjArr + +-- | Get all vertices of the graph. +-- +getVertices :: GenericGraph v e -> [v] +getVertices (GenericGraph idxArr _ _) = map snd $ A.assocs idxArr + +-- | Get subgraph on given vertices. Note that indexation will be CHANGED. +-- Be careful with !. and ?. operators. +-- +subgraph :: Ord v => GenericGraph v e -> [Int] -> GenericGraph v e +subgraph graph toKeep = fromList (newVertices, newEdges) + where + vSet :: S.Set Int + vSet = S.fromList toKeep + + eRemain :: (Int, Int, e) -> Bool + eRemain (at, other, _) = (at `S.member` vSet) && (other `S.member` vSet) + + (oldVertices, edges) = filter eRemain <$> toList graph + (newVertices, oldIdx) = unzip . filter (\(_, ix) -> ix `S.member` vSet) $ zip oldVertices [0..] + + vMap :: Map Int Int + vMap = M.fromList $ zip oldIdx [0 ..] + + newEdges = map (\(at, other, bond) -> (vMap ! at, vMap ! other, bond)) edges + +-- | Add given vertices to graph. +-- +addVertices :: Ord v => GenericGraph v e -> [v] -> GenericGraph v e +addVertices graph toAdd = fromList (first (++ toAdd) (toList graph)) + +-- | Remove given vertices from the graph. Note that indexation will be CHANGED. +-- Be careful with !. and ?. operators. +-- +removeVertices :: Ord v => GenericGraph v e -> [Int] -> GenericGraph v e +removeVertices graph toRemove = fromList (newVertices, newEdges) + where + vSet :: S.Set Int + vSet = S.fromList toRemove + + eRemove :: (Int, Int, e) -> Bool + eRemove (at, other, _) = (at `S.notMember` vSet) && (other `S.notMember` vSet) + + (oldVertices, edges) = filter eRemove <$> toList graph + (newVertices, oldIdx) = unzip . filter ((`S.notMember` vSet) . snd) $ zip oldVertices [0..] + + vMap :: Map Int Int + vMap = M.fromList $ zip oldIdx [0 ..] + + newEdges = map (\(at, other, bond) -> (vMap ! at, vMap ! other, bond)) edges + +-- | Remove given edges from the graph. Note that isolated vertices are allowed. +-- This will NOT affect indexation. +-- +removeEdges :: Ord v => GenericGraph v e -> [(Int, Int)] -> GenericGraph v e +removeEdges graph toRemove = fromList (vertices, edges) + where + eSet :: S.Set (Int, Int) + eSet = S.fromList toRemove + + (vertices, edges) = filter eRemove <$> toList graph + + eRemove (at, other, _) = ((at, other) `S.notMember` eSet) && ((other, at) `S.notMember` eSet) + +-- | Add given edges to the graph. +-- +addEdges :: Ord v => GenericGraph v e -> [(Int, Int, e)] -> GenericGraph v e +addEdges (GenericGraph inds rinds edges) toAdd = GenericGraph inds rinds res + where + accumList = foldl (\x (a, b, t) -> x ++ [(a, (b, t)), (b, (a, t))]) [] toAdd + res = A.accum (flip (:)) edges accumList + +-- | Returns type of edge with given starting and ending indices. +-- +typeOfEdge :: Ord v => GenericGraph v e -> Int -> Int -> e +typeOfEdge graph fromInd toInd = res + where + neighbors = gAdjacency graph A.! fromInd + res = snd (fromJust (find ((== toInd) . fst) neighbors)) + +-- | Safe extraction from the graph. If there is no requested key in it, +-- empty list is returned. +-- +safeIdx :: GenericGraph v e -> Int -> [Int] +safeIdx graph = map fst . fromMaybe [] . (graph ?.) + +-- | Safe extraction from the graph. If there is no requested key in it, +-- empty list is returned. +-- +safeAt :: GenericGraph v e -> Int -> [(Int, e)] +safeAt graph = fromMaybe [] . (graph ?.) + +-- | Get edge from graph, which starting and ending indices match +-- given indices. +-- +getEdge :: GenericGraph v e -> Int -> Int -> e +getEdge graph from to = found + where + neighbors = graph !. from + found = snd (fromJust (find ((== to) . fst) neighbors)) + +-- | Check that two vertices with given indexes have edge between them. +-- +isConnected :: GenericGraph v e -> Int -> Int -> Bool +isConnected g fInd tInd = isJust $ find ((==) tInd . fst) $ safeAt g fInd + +-- | Returns graph that is the sum of two given graphs assuming that they are disjoint. +-- +sumGraphs :: Ord v => GenericGraph v e -> GenericGraph v e -> GenericGraph v e +sumGraphs graphA graphB = res + where + (vertA, edgeA) = toList graphA + (vertB, edgeB) = toList graphB + renameMapB = M.fromList (zip [0..length vertB - 1] [length vertA..length vertA + length vertB - 1]) + renameFunc = (renameMapB M.!) + + newVertices = vertA ++ vertB + newEdges = edgeA ++ fmap (\(a, b, t) -> (renameFunc a, renameFunc b, t)) edgeB + + res = fromList (newVertices, newEdges) diff --git a/src/Math/Grads/Graph.hs b/src/Math/Grads/Graph.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Graph.hs @@ -0,0 +1,92 @@ +-- | Module that provides 'Graph' type class and several useful functions +-- for interaction with 'Graph's. +-- +module Math.Grads.Graph + ( EdgeList + , Graph (..) + , GraphEdge + , changeIndsEdge + , changeTypeEdge + , edgeType + ) where + +import Data.List (nub) + +-- | 'GraphEdge' is just triple, containing index of starting vertex of edge, +-- index of ending vertex of edge and edge's type. +-- +type GraphEdge e = (Int, Int, e) + +-- | Type alias for list of 'GraphEdge's. +-- +type EdgeList e = [GraphEdge e] + +-- | Get edge's type from 'GraphEdge'. +-- +edgeType :: GraphEdge e -> e +edgeType (_, _, t) = t + +-- | Given transformation of edge types transforms 'GraphEdge'. +-- +changeTypeEdge :: (e1 -> e2) -> GraphEdge e1 -> GraphEdge e2 +changeTypeEdge f (a, b, t) = (a, b, f t) + +-- | Given transformation of edge's indices transforms 'GraphEdge'. +-- +changeIndsEdge :: (Int -> Int) -> GraphEdge e -> GraphEdge e +changeIndsEdge f (a, b, t) = (f a, f b, t) + +-- | Type class that gives data structure properties of graph. +-- +class Graph g where + -- | Construct a graph from list of vertices and edges. + -- + fromList :: (Ord v, Eq v) => ([v], [GraphEdge e]) -> g v e + + -- | Get a list of all vertices and edges from the graph. + -- + toList :: (Ord v, Eq v) => g v e -> ([v], [GraphEdge e]) + + -- | Get the number of vertices. + -- + vCount :: g v e -> Int + + -- | Unsafe get adjacent vertices. + -- + infixl 9 !> + (!>) :: (Ord v, Eq v) => g v e -> v -> [(v, e)] + + -- | Unsafe get adjacent indices. + -- + infixl 9 !. + (!.) :: g v e -> Int -> [(Int, e)] + + -- | Safe get adjacent vertices. + -- + infixl 9 ?> + (?>) :: (Ord v, Eq v) => g v e -> v -> Maybe [(v, e)] + + -- | Safe get adjacent indices. + -- + infixl 9 ?. + (?.) :: g v e -> Int -> Maybe [(Int, e)] + + -- | Get a list of edges starting at given vertex. + -- + incident :: (Ord v, Eq v) => g v e -> v -> [(v, v, e)] + incident gr at = (\(a, b) -> (at, a, b)) <$> gr !> at + + -- | Safe get a list of edges starting at given vertex. + -- + safeIncident :: (Ord v, Eq v) => g v e -> v -> Maybe [(v, v, e)] + safeIncident gr at = map (\(a, b) -> (at, a, b)) <$> gr ?> at + + -- | Get a list of index edges starting at given vertex. + -- + incidentIdx :: (Eq e) => g v e -> Int -> [GraphEdge e] + incidentIdx gr idx = nub ((\(a, b) -> (min idx a, max idx a, b)) <$> gr !. idx) + + -- | Safe get a list of index edges starting at given vertex. + -- + safeIncidentIdx :: (Eq e) => g v e -> Int -> Maybe [GraphEdge e] + safeIncidentIdx gr idx = nub <$> (map (\(a, b) -> (min idx a, max idx a, b)) <$> gr ?. idx) diff --git a/src/Math/Grads/Utils.hs b/src/Math/Grads/Utils.hs new file mode 100644 --- /dev/null +++ b/src/Math/Grads/Utils.hs @@ -0,0 +1,26 @@ +-- | Different utility functions for usage in Math.Grads. +-- +module Math.Grads.Utils + ( nub + , subsets + , uniter + ) where + +import Data.List (group, sort) + +-- | nub that works in O(n log n) time. +-- +nub :: (Ord a, Eq a) => [a] -> [a] +nub = fmap head . group . sort + +-- | Zips list with its tail. +-- +uniter :: [a] -> [(a, a)] +uniter [] = [] +uniter l = zip l $ drop 1 l + +-- | Returns all possible subsets of given list as list of lists. +-- +subsets :: [a] -> [[a]] +subsets [] = [[]] +subsets (x : xs) = subsets xs ++ ((x :) <$> subsets xs) diff --git a/test/Coords.hs b/test/Coords.hs new file mode 100644 --- /dev/null +++ b/test/Coords.hs @@ -0,0 +1,63 @@ +{-# LANGUAGE MultiParamTypeClasses #-} +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Main where + +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M +import Math.Grads.Drawing.Coords (Drawable, getCoordsForGraph) +import Math.Grads.GenericGraph (GenericGraph) +import Math.Grads.Graph (fromList) +import System.Random (mkStdGen) +import Test.Hspec + +instance Drawable GenericGraph Int Int + +pathToGraphs :: FilePath +pathToGraphs = "data/Graphs.txt" + +roundPair :: (Float, Float) -> (Int, Int) +roundPair (a, b) = (round a, round b) + +testMap :: IO (Map String (GenericGraph Int Int, Map Int (Int, Int))) +testMap = do + graphsInLines <- lines <$> readFile pathToGraphs + let graphsInWords = fmap words graphsInLines + + let forMap = fmap (\(x : y : z : _) -> (x, (fromList (read y), fmap roundPair (read z)))) graphsInWords + return (M.fromList forMap) + +testDrawing :: SpecWith () +testDrawing = describe "Check whether molecules are being drawn correctly." $ do + it "Only path" $ do + (graph, coords) <- fmap (M.! "only_path") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Just coords + it "Only cycles" $ do + (graph, coords) <- fmap (M.! "only_cycles") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Just coords + it "Simple drawing" $ do + (graph, coords) <- fmap (M.! "simple_drawing") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Just coords + it "Hard drawing" $ do + (graph, coords) <- fmap (M.! "hard_drawing") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Just coords + it "Paths through conjugated cycles" $ do + (graph, coords) <- fmap (M.! "paths_through_conjugated_cycles") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Just coords + +testErrors :: SpecWith () +testErrors = describe "Check that coordinates for molecules that we can't draw are returned as Nothing." $ do + it "Too big cycle" $ do + (graph, _) <- fmap (M.! "too_big_cycle") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Nothing + it "Bad conjugated cycle" $ do + (graph, _) <- fmap (M.! "bad_conjugated_cycle") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Nothing + it "Disappearing cycle" $ do + (graph, _) <- fmap (M.! "disappearing_cycle") testMap + (roundPair <$>) <$> getCoordsForGraph (mkStdGen 0) graph `shouldBe` Nothing + +main :: IO () +main = hspec $ do + testDrawing + testErrors diff --git a/test/Graph.hs b/test/Graph.hs new file mode 100644 --- /dev/null +++ b/test/Graph.hs @@ -0,0 +1,54 @@ +module Main where + +import Data.List (sort) +import Math.Grads.GenericGraph (GenericGraph, removeEdges, + removeVertices, subgraph) +import Math.Grads.Graph (fromList, incident, safeIncident, + (!.), (!>), (?>)) +import Test.Hspec + +main :: IO () +main = hspec $ do + opTests + subgraphTests + removalTests + +vertices :: [Int] +vertices = [0, 1, 2, 5, 7, 8, 9, 10, 11, 19, 31, 78, 79, 85, 99, 60, 53, 100, 42] + +edges :: [(Int, Int, Int)] +edges = [(0, 1, 1), (1, 2, 3), (1, 1, 2), (0, 2, 4), (3, 4, 2), (3, 5, 1), + (5, 6, 2), (5, 7, 5), (7, 8, 3), (0, 9, 7), (0, 10, 2), (10, 11, 2), + (10, 13, 565), (13, 14, 546), (15, 14, 42), (15, 16, -4), (16, 17, 0), + (18, 16, 1), (16, 12, 0), (12, 10, 1)] + +graph :: GenericGraph Int Int +graph = fromList (vertices, edges) + +opTests :: Spec +opTests = describe "Operations on graphs." $ do + it "Adjacent to 0." $ sort (graph !> 0) `shouldBe` [(1, 1), (2, 4), (19, 7), (31, 2)] + it "Adjacent to 1." $ sort <$> (graph ?> 1) `shouldBe` Just [(0, 1), (1, 2), (2, 3)] + it "Adjacent to 5." $ sort (graph !> 5) `shouldBe` [(7, 2), (8, 1)] + it "Adjacent to 14." $ (graph ?> 14) `shouldBe` Nothing + it "Edges incident to 8." $ sort (graph `incident` 8) `shouldBe` [(8, 5, 1), (8, 9, 2), (8, 10, 5)] + it "Edges incident to 53." $ sort <$> (graph `safeIncident` 53) `shouldBe` Just [(53, 42, 1), (53, 60, -4), (53, 79, 0), (53, 100, 0)] + +subgraphTests :: Spec +subgraphTests = describe "Subgraph tests." $ do + let subg = graph `subgraph` [0, 3, 5, 7, 8, 11, 14] + it "Adjacent to 0." $ subg !> 0 `shouldBe` [] + it "Adjacent to 3." $ subg !. 1 `shouldBe` [(2, 1)] + it "Adjacent to 8." $ sort (subg !. 2) `shouldBe` [(1, 1), (3, 5)] + +removalTests :: Spec +removalTests = describe "Remove operations tests." $ do + let g1 = graph `removeEdges` [(5, 6), (15, 14), (0, 10), (10, 12)] + let g2 = graph `removeVertices` [1, 3, 5, 7, 10, 15] + it "Adjacent to 0." $ sort (g1 !> 0) `shouldBe` [(1, 1), (2, 4), (19, 7)] + it "Adjacent to 31." $ sort (g1 !> 31) `shouldBe` [(78, 2), (85, 565)] + it "Edges incident to 8." $ sort (g1 `incident` 8) `shouldBe` [(8, 5, 1), (8, 10, 5)] + it "Adjacent to 14." $ (g1 ?> 14) `shouldBe` Nothing + it "Adjacent to 0." $ sort (g2 !> 0) `shouldBe` [(2, 4), (19, 7)] + it "Adjacent to 10." $ g2 ?> 31 `shouldBe` Nothing + it "Adjacent to 11." $ g2 !> 78 `shouldBe` [] diff --git a/test/Isomorphism.hs b/test/Isomorphism.hs new file mode 100644 --- /dev/null +++ b/test/Isomorphism.hs @@ -0,0 +1,96 @@ +{-# OPTIONS_GHC -fno-warn-orphans #-} +{-# LANGUAGE MultiParamTypeClasses #-} + +module Main where + +import qualified Data.Array as A +import Data.Map.Strict (Map) +import qualified Data.Map.Strict as M +import Math.Grads.Algo.Isomorphism (GComparable (..), isIsoSub) +import Math.Grads.GenericGraph (GenericGraph, gIndex) +import Math.Grads.Graph (fromList) +import Test.Hspec + +instance GComparable GenericGraph Int Int GenericGraph Int Int where + vComparator g1 g2 ind1 ind2 = gIndex g1 A.! ind1 == gIndex g2 A.! ind2 + eComparator _ _ (_, _, t) (_, _, t') = t == t' + +pathToGraphs :: FilePath +pathToGraphs = "data/Graphs.txt" + +testMap :: IO (Map String (GenericGraph Int Int)) +testMap = do + graphsInLines <- lines <$> readFile pathToGraphs + let graphsInWords = fmap words graphsInLines + + let forMap = fmap (\(x : y : _) -> (x, fromList (read y))) graphsInWords + return (M.fromList forMap) + +bigSubGraph :: GenericGraph Int Int +bigSubGraph = fromList ( [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + , [ (0, 1, 1), (0, 26, 1), (1, 2, 1), (2, 3, 1), (3, 4, 1), (3, 25, 1), (4, 5, 1), (4, 11, 1), (5, 6, 1) + , (6, 7, 1), (7, 8, 1), (7, 10, 1), (8, 9, 1), (9, 10, 1), (11, 12, 1), (11, 24, 1), (12, 13, 1), (13, 14, 1) + , (13, 18, 1), (14, 15, 1), (15, 16, 1), (16, 17, 1), (17, 18, 1), (18, 19, 1), (19, 20, 1), (19, 24, 1) + , (20, 21, 1), (21, 22, 1), (22, 23, 1), (23, 24, 1), (25, 26, 1) + ] + ) + +pathGraph :: GenericGraph Int Int +pathGraph = fromList ([0, 0, 0, 0, 0, 0, 0], [(0, 1, 1), (0, 2, 1), (0, 3, 1), (0, 4, 1), (4, 5, 1), (4, 6, 1)]) + +conjugatedCycles :: GenericGraph Int Int +conjugatedCycles = fromList ( [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + , [ (0, 1, 1), (0, 2, 1), (1, 3, 1), (2, 4, 1), (4, 5, 1), (3, 5, 1), (3, 6, 1), (5, 7, 1) + , (6, 8, 1), (7, 9, 1), (8, 9, 1), (1, 10, 1), (6, 11, 1), (10, 12, 1), (11, 12, 1) + ] + ) + +connectedCycles :: GenericGraph Int Int +connectedCycles = fromList ( [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + , [ (0, 1, 1), (0, 2, 1), (1, 3, 1), (2, 4, 1), (4, 5, 1), (3, 5, 1), (3, 6, 1), (6, 7, 1) + , (6, 8, 1), (7, 9, 1), (8, 10, 1), (9, 11, 1), (10, 11, 1), (8, 12, 1) + ] + ) + +cycleAndTriangle :: GenericGraph Int Int +cycleAndTriangle = fromList ( [0, 0, 0, 0, 0, 0, 0, 0, 0, 0] + , [ (0, 1, 1), (0, 2, 1), (1, 3, 1), (2, 4, 1), (3, 4, 1), (4, 5, 1), (5, 6, 1), (6, 7, 1) + , (7, 8, 1), (7, 9, 1), (8, 9, 1) + ] + ) + +triangleAndTriangle :: GenericGraph Int Int +triangleAndTriangle = fromList ( [0, 0, 0, 0, 0, 0, 0] + , [(0, 1, 1), (0, 2, 1), (1, 2, 1), (1, 3, 1), (3, 4, 1), (3, 5, 1), (4, 5, 1)] + ) + +testIsIsoSub :: SpecWith () +testIsIsoSub = describe "Check whether subgraph isomorphism algorithm is working correctly" $ do + it "Path" $ do + graph <- fmap (M.! "only_path") testMap + graph `shouldSatisfy` isIsoSub pathGraph + it "Conjugated cycles" $ do + graph <- fmap (M.! "only_cycles") testMap + graph `shouldSatisfy` isIsoSub conjugatedCycles + it "Connected cycles" $ do + graph <- fmap (M.! "simple_drawing") testMap + graph `shouldSatisfy` isIsoSub connectedCycles + it "Conjugated cycles again" $ do + graph <- fmap (M.! "hard_drawing") testMap + graph `shouldSatisfy` isIsoSub conjugatedCycles + it "Cycle and triangle" $ do + graph <- fmap (M.! "paths_through_conjugated_cycles") testMap + graph `shouldSatisfy` isIsoSub cycleAndTriangle + it "Big graph" $ do + graph <- fmap (M.! "takes_long_if_done_wrong") testMap + graph `shouldSatisfy` isIsoSub bigSubGraph + it "Triangle and triangle. No match" $ do + graph <- fmap (M.! "paths_through_conjugated_cycles") testMap + graph `shouldNotSatisfy` isIsoSub triangleAndTriangle + it "Cycle and triangle. No match" $ do + graph <- fmap (M.! "simple_drawing") testMap + graph `shouldNotSatisfy` isIsoSub cycleAndTriangle + +main :: IO () +main = hspec $ do + testIsIsoSub