RankingTree.h

Go to the documentation of this file.

#ifndef INCLUDE_MOLASSEMBLER_RANKING_HIERARCHICAL_TREE_H
#define INCLUDE_MOLASSEMBLER_RANKING_HIERARCHICAL_TREE_H

#include "boost/variant.hpp"

#include "Molassembler/Detail/BuildTypeSwitch.h"
#include "Molassembler/AtomStereopermutator.h"
#include "Molassembler/Modeling/BondDistance.h"
#include "Molassembler/BondStereopermutator.h"
#include "Molassembler/Log.h"
#include "Molassembler/Molecule.h"
#include "Molassembler/Molecule/OrderDiscoveryHelper.h"

#include "Molassembler/Temple/Adaptors/AllPairs.h"

using namespace std::string_literals;

namespace Scine {
namespace Molassembler {

class RankingTree {
private:
  static unsigned debugMessageCounter_;

  static std::string adaptMolGraph_(std::string molGraph);

  static void writeGraphvizFiles_(
    const std::vector<std::string>& graphvizStrings
  );

public:

  enum class ExpansionOption {
    OnlyRequiredBranches,
    Full
  };

  struct VertexData {
    AtomIndex molIndex;
    bool isDuplicate;

    boost::optional<AtomStereopermutator> stereopermutatorOption;

    // This is the place for eventual atomic number deviations
    // boost::optional<double> deviantAtomicNumber;
  };

  struct EdgeData {
    boost::optional<BondStereopermutator> stereopermutatorOption;
  };

  using BglType = boost::adjacency_list<
    /* OutEdgeListS = Type of Container for edges of a vertex
     * Options: vector, list, slist, set, multiset, unordered_set
     * Choice: setS
     */
    boost::setS,
    /* VertexListS = Type of Container for vertices
     * Options: vector, list, slist, set, multiset, unordered_set
     * Choice: vecS, removing vertices does not occur
     */
    boost::vecS,
    /* DirectedS = Is the graph directed or not?
     * Choice: Bidirectional. Although a tree is merely directional, and not
     *   bi-directional, we nevertheless also need fast access to in_edges, so
     *   we choose bidirectional
     */
    boost::bidirectionalS,
    // VertexProperty = What information is stored about vertices?
    VertexData,
    // EdgeProperty = What information is stored about edges?
    EdgeData
  >;

  // IUPAC Sequence rule one tree vertex comparator
  class SequenceRuleOneVertexComparator;

  // IUPAC Sequence rule two tree vertex comparator
  class SequenceRuleTwoVertexComparator;

  // IUPAC Sequence rule three tree edge comparator
  class SequenceRuleThreeEdgeComparator;

  // IUPAC Sequence rule four tree vertex and edge mixed comparator
  class SequenceRuleFourVariantComparator;

  // IUPAC Sequence rule five tree vertex and edge mixed comparator
  class SequenceRuleFiveVariantComparator;

  // Returns whether the vertex or edge has an instantiated Stereopermutator
  class VariantHasInstantiatedStereopermutator;

  // Returns the mixed depth (see mixedDepth_ functions) of the vertex or edge
  class VariantDepth;

  // Returns the source node of a vertex (identity) or an edge (source node)
  class VariantSourceNode;

  // Returns a string representation of the stereopermutator on vertex or edge
  class VariantStereopermutatorStringRepresentation;

  /*
   * Returns whether two variants' stereopermutators form a like pair or not
   * according to the sub-rules in IUPAC sequence rule four.
   */
  class VariantLikePair;

  class GraphvizWriter;

private:
  using TreeVertexIndex = BglType::vertex_descriptor;
  using TreeEdgeIndex = BglType::edge_descriptor;
  using VariantType = boost::variant<TreeVertexIndex, TreeEdgeIndex>;

  struct JunctionInfo {
    TreeVertexIndex junction;

    std::vector<TreeVertexIndex> firstPath, secondPath;
  };

  static constexpr TreeVertexIndex rootIndex = 0;

/* State */
  BglType tree_;

  OrderDiscoveryHelper<TreeVertexIndex> branchOrderingHelper_;

  // Closures
  const Graph& graph_;
  const StereopermutatorList& stereopermutatorsRef_;
  const std::string adaptedMolGraphviz_;

/* Minor helper classes and functions */
  TreeVertexIndex parent_(const TreeVertexIndex& index) const;

  std::vector<TreeVertexIndex> adjacents_(TreeVertexIndex index) const;

  unsigned adjacentTerminalHydrogens_(const TreeVertexIndex& index) const;

  bool isBondSplitDuplicateVertex_(const TreeVertexIndex& index) const;

  bool isCycleClosureDuplicateVertex_(const TreeVertexIndex& index) const;

  std::vector<TreeVertexIndex> auxiliaryAdjacentsToRank_(
    TreeVertexIndex sourceIndex,
    const std::vector<TreeVertexIndex>& excludeIndices
  ) const;

  template<typename ComparisonSets>
  bool notEmpty_(const ComparisonSets& comparisonSets) const {
    return Temple::all_of(
      comparisonSets,
      [](const auto& mapPair) -> bool {
        return !mapPair.second.empty();
      }
    );
  }

  unsigned nonDuplicateDegree_(const TreeVertexIndex& index) const;

  std::vector<TreeVertexIndex> addBondOrderDuplicates_(
    const TreeVertexIndex& treeSource,
    const TreeVertexIndex& treeTarget
  );

  std::unordered_set<TreeVertexIndex> treeIndicesInBranch_(TreeVertexIndex index) const;

  std::unordered_set<AtomIndex> molIndicesInBranch_(TreeVertexIndex index) const;

  unsigned duplicateDepth_(TreeVertexIndex index) const;

  unsigned depthOfNode_(TreeVertexIndex index) const;

  unsigned mixedDepth_(const TreeVertexIndex& vertexIndex) const;

  unsigned mixedDepth_(const TreeEdgeIndex& edgeIndex) const;

  JunctionInfo junction_(const TreeVertexIndex& a, const TreeVertexIndex& b) const;

  bool molIndexExistsInBranch_(
    AtomIndex molIndex,
    TreeVertexIndex treeIndex
  ) const;

/* Major helper functions */
  std::vector<TreeVertexIndex> expand_(
    const TreeVertexIndex& index,
    const std::unordered_set<AtomIndex>& molIndicesInBranch
  );

  template<typename SetValueType, typename ComparatorType>
  bool multisetCompare_(
    const std::multiset<SetValueType, ComparatorType>& a,
    const std::multiset<SetValueType, ComparatorType>& b
  ) const {
    return std::lexicographical_compare(
      std::begin(a),
      std::end(a),
      std::begin(b),
      std::end(b),
      ComparatorType {*this}
    );
  }

  template<typename SetValueType, typename ComparatorType>
  void compareBFSSets_(
    const std::map<
      TreeVertexIndex,
      std::multiset<SetValueType, ComparatorType>
    >& comparisonSets,
    const std::vector<
      std::vector<TreeVertexIndex>
    >& undecidedSets,
    OrderDiscoveryHelper<TreeVertexIndex>& orderingHelper
  ) const {
    for(const auto& undecidedSet : undecidedSets) {
      Temple::forEach(
        Temple::Adaptors::allPairs(undecidedSet),
        [&](const TreeVertexIndex a, const TreeVertexIndex b) {
          if(multisetCompare_(comparisonSets.at(a), comparisonSets.at(b))) {
            orderingHelper.addLessThanRelationship(a, b);
          } else if(multisetCompare_(comparisonSets.at(b), comparisonSets.at(a))) {
            orderingHelper.addLessThanRelationship(b, a);
          }
        }
      );
    }
  }

  static std::string toString(TreeVertexIndex vertex);
  std::string toString(const TreeEdgeIndex& edge) const;
  std::string toString(const VariantType& vertexOrEdge) const;

  template<
    unsigned ruleNumber,
    bool bfsDownOnly,
    bool insertEdges,
    bool insertVertices,
    typename MultisetValueType,
    typename MultisetComparatorType
  > void runBFS_(
    TreeVertexIndex sourceIndex,
    OrderDiscoveryHelper<TreeVertexIndex>& orderingHelper,
    const boost::optional<unsigned>& depthLimitOptional = boost::none
  ) const {
    /* Static safety checks */
    static_assert(
      insertEdges || insertVertices,
      "During BFS, you must insert either edges, vertices, or both into the "
      " comparison sets, but definitely not neither"
    );

    static_assert(
      (
        insertEdges
        && !std::is_same<MultisetValueType, TreeVertexIndex>::value
      ) || (
        insertVertices
        && !std::is_same<MultisetValueType, TreeEdgeIndex>::value
      ),
      "Multiset value type and insert booleans are mismatched"
    );

    // Check if depthLimit prohibits first comparison too
    if(depthLimitOptional.value_or(std::numeric_limits<unsigned>::max()) == 0) {
      return;
    }

    /* Declarations */

    /* For each branch to compare, keep a multiset of all the encountered
     * indices in the branch. The multiset keeps a descending order of the
     * indices according to sequence rule 1 and can be lexicographically
     * compared against one another using the SequenceRuleOneVertexComparator
     *
     * NOTE: Do not use operator [] to access members in comparisonSets. This
     * old form of accessible and assignable proxy object forces a
     * default-instantiation of the Comparator, which the
     * MultisetComparatorTypes cannot do. Always use .at().
     */
    std::map<
      TreeVertexIndex,
      std::multiset<MultisetValueType, MultisetComparatorType>
    > comparisonSets;

    /* For each branch to compare, keep a set of seed indices to follow in an
     * iteration. These are expanded in each iteration as long as the branch
     * they contain remains relevant (i.e. its relation to the other branches
     * is still unclear): They themselves are placed in the comparisonSet,
     * while their respective adjacents are the new seeds for the next
     * iteration.
     */
    std::map<
      TreeVertexIndex,
      std::vector<TreeVertexIndex>
    > seeds;

    /* In case the BFS in not down-only, we have to track which indices we
     * have visited to ensure BFS doesn't backtrack or reuse vertices.
     */
    std::unordered_set<TreeVertexIndex> visitedVertices;


    /* Initialization */
    if(!bfsDownOnly) { // Initialization is only necessary in this case
      visitedVertices.insert(sourceIndex);
    }

    auto undecidedSets = orderingHelper.getUndecidedSets();

    for(const auto& undecidedSet : undecidedSets) {
      for(const auto& undecidedBranch : undecidedSet) {
        comparisonSets.emplace(
          undecidedBranch,
          *this
        );

        if constexpr(insertVertices) {
          comparisonSets.at(undecidedBranch).insert(undecidedBranch);
        }

        if constexpr(insertEdges) {
          if(bfsDownOnly) {
            const auto edge = boost::edge(sourceIndex, undecidedBranch, tree_).first;
            comparisonSets.at(undecidedBranch).insert(edge);
          } else {
            // Need to be direction-agnostic
            auto forwardEdge = boost::edge(sourceIndex, undecidedBranch, tree_);
            auto backwardEdge = boost::edge(undecidedBranch, sourceIndex, tree_);

            const auto edge = forwardEdge.second ? forwardEdge.first : backwardEdge.first;
            comparisonSets.at(undecidedBranch).insert(edge);
          }
        }

        seeds[undecidedBranch] = {undecidedBranch};
      }
    }

    // Compare undecided multisets, and add any discoveries to the ordering helper
    compareBFSSets_(comparisonSets, undecidedSets, orderingHelper);

    // Update the undecided sets
    undecidedSets = orderingHelper.getUndecidedSets();

    if constexpr (buildTypeIsDebug) {
      // Write debug graph files if the corresponding log particular is set
      if(
        Log::particulars.count(Log::Particulars::RankingTreeDebugInfo) > 0
        && notEmpty_(comparisonSets)
      ) {
        std::string header = (
          (
            bfsDownOnly
            ? ""s
            : "aux "s
          ) + "R"s + std::to_string(ruleNumber)
        );

        writeGraphvizFiles_({
          adaptedMolGraphviz_,
          dumpGraphviz(
            header,
            {sourceIndex},
            collectSeeds_(seeds, undecidedSets)
          ),
          makeBFSStateGraph_(
            header,
            sourceIndex,
            comparisonSets,
            undecidedSets
          ),
          orderingHelper.dumpGraphviz()
        });
      }
    }

    unsigned depth = 1;

    /* Loop BFS */

    /* As long as there are undecided sets and seeds whose expansion could be
     * relevant for those undecided sets, continue BFS
     */
    while(
      // Undecided indices remain
      !undecidedSets.empty()
      // Seeds exist that are relevant to the undecided sets
      && relevantSeeds_(seeds, undecidedSets)
      && depth < depthLimitOptional.value_or(std::numeric_limits<unsigned>::max())
    ) {
      // BFS Step
      for(const auto& undecidedSet: undecidedSets) {
        for(const auto& undecidedBranch: undecidedSet) {
          /* Clear the comparison set at undecided branches prior to inserting
           * anything. This doesn't change anything semantically, since
           * everything from the prior iteration (or initialization) must have
           * compared equal within this undecidedSet.
           *
           * Makes for cleaner visualization and reduces the number of
           * Comparator calls during insert and lexicographical_compare.
           */
          comparisonSets.at(undecidedBranch).clear();

          std::vector<TreeVertexIndex> newSeeds;

          for(const auto& seed : seeds.at(undecidedBranch)) {

            if(!bfsDownOnly) {
              // Mark as visited
              visitedVertices.insert(seed);

              // In-edges are only relevant for omnidirectional BFS
              for(
                auto inIterPair = boost::in_edges(seed, tree_);
                inIterPair.first != inIterPair.second;
                ++inIterPair.first
              ) {
                const auto& inEdge = *inIterPair.first;

                auto edgeSource = boost::source(inEdge, tree_);

                // Check if already placed
                if(visitedVertices.count(edgeSource) == 0) {
                  if constexpr(insertVertices) {
                    comparisonSets.at(undecidedBranch).insert(edgeSource);
                  }

                  if constexpr(insertEdges) {
                    comparisonSets.at(undecidedBranch).insert(inEdge);
                  }

                  newSeeds.push_back(edgeSource);
                }
              }
            }

            // Out edges are considered in all cases
            for( // Out-edges
              auto outIterPair = boost::out_edges(seed, tree_);
              outIterPair.first != outIterPair.second;
              ++outIterPair.first
            ) {
              const auto& outEdge = *outIterPair.first;

              auto edgeTarget = boost::target(outEdge, tree_);

              // Skip this vertex if in omnidirectional BFS and already-seen node
              if(!bfsDownOnly && visitedVertices.count(edgeTarget) > 0) {
                continue;
              }

              if constexpr(insertVertices) {
                comparisonSets.at(undecidedBranch).insert(edgeTarget);
              }

              if constexpr(insertEdges) {
                comparisonSets.at(undecidedBranch).insert(outEdge);
              }

              // Add out edge target to seeds only if non-terminal
              if(!tree_[edgeTarget].isDuplicate) {
                newSeeds.push_back(edgeTarget);
              }
            }
          }

          // Overwrite seeds
          seeds.at(undecidedBranch) = std::move(newSeeds);
        }
      }

      // Make comparisons in all undecided sets
      compareBFSSets_(comparisonSets, undecidedSets, orderingHelper);

      // Recalculate the undecided sets
      undecidedSets = orderingHelper.getUndecidedSets();

      // Increment depth
      ++depth;

      if constexpr (buildTypeIsDebug) {
        if(
          Log::particulars.count(Log::Particulars::RankingTreeDebugInfo) > 0
          && notEmpty_(comparisonSets)
        ) {
          std::string header;
          if(!bfsDownOnly) {
            header += "aux "s;
          }
          header += "R"s + std::to_string(ruleNumber);

          writeGraphvizFiles_({
            adaptedMolGraphviz_,
            dumpGraphviz(
              header,
              {sourceIndex},
              collectSeeds_(seeds, undecidedSets)
            ),
            makeBFSStateGraph_(
              header,
              sourceIndex,
              comparisonSets,
              undecidedSets
            ),
            orderingHelper.dumpGraphviz()
          });
        }
      }
    }
  }

  // Flatten the undecided branches' seeds into a single set
  static std::unordered_set<TreeVertexIndex> collectSeeds_(
    const std::map<
      TreeVertexIndex,
      std::vector<TreeVertexIndex>
    >& seeds,
    const std::vector<
      std::vector<TreeVertexIndex>
    >& undecidedSets
  );

  template<typename ValueType, typename ComparatorType>
  std::string makeBFSStateGraph_(
    const std::string& title,
    const TreeVertexIndex& base,
    const std::map<
      TreeVertexIndex,
      std::multiset<ValueType, ComparatorType>
    >& comparisonSet,
    const std::vector<
      std::vector<TreeVertexIndex>
    >& undecidedSets
  ) const {

    std::set<TreeVertexIndex> activeBranches;
    for(const auto& set : undecidedSets) {
      for(const auto& value : set) {
        activeBranches.insert(value);
      }
    }

    struct DisplayGraphVertexData {
      std::string representation;
    };

    using DisplayGraph = boost::adjacency_list<
      /* OutEdgeListS = Type of Container for edges of a vertex
       * Options: vector, list, slist, set, multiset, unordered_set
       * Choice: setS
       */
      boost::setS,
      /* VertexListS = Type of Container for vertices
       * Options: vector, list, slist, set, multiset, unordered_set
       * Choice: vecS, removing vertices does not occur
       */
      boost::vecS,
      /* DirectedS = Is the graph directed or not?
       * Choice: Directed, do not need bidirectional
       */
      boost::directedS,
      /* VertexProperty = What information is stored about vertices?
       * Choice: Atom, containing an index and an element type
       */
      DisplayGraphVertexData
    >;

    std::set<TreeVertexIndex> colorVertices;
    std::set<TreeVertexIndex> squareVertices;

    DisplayGraph graph;

    auto baseVertex = boost::add_vertex(graph);
    graph[baseVertex].representation = toString(base);

    squareVertices.insert(baseVertex);

    for(const auto& branchIterPair : comparisonSet) {
      const auto& treeBranchIndex = branchIterPair.first;
      const auto& treeBranchMultiset = branchIterPair.second;

      auto branchVertex = boost::add_vertex(graph);
      graph[branchVertex].representation = toString(treeBranchIndex);

      squareVertices.insert(branchVertex);

      if(activeBranches.count(treeBranchIndex) > 0) {
        colorVertices.insert(branchVertex);
      }

      boost::add_edge(baseVertex, branchVertex, graph);

      typename DisplayGraph::vertex_descriptor continuation = branchVertex;
      for(const auto& multisetValue : treeBranchMultiset) {
        auto newVertex = boost::add_vertex(graph);
        graph[newVertex].representation = toString(multisetValue);

        boost::add_edge(continuation, newVertex, graph);
        continuation = newVertex;
      }
    }

    class SmallGraphWriter {
    private:
      const DisplayGraph& graphBase_;
      const std::string title_;
      const std::set<TreeVertexIndex> colorVertices_;
      const std::set<TreeVertexIndex> squareVertices_;

    public:
      SmallGraphWriter(
        const DisplayGraph& base,
        std::string title,
        std::set<TreeVertexIndex> colorVertices,
        std::set<TreeVertexIndex> squareVertices
      ) : graphBase_(base),
          title_(std::move(title)),
          colorVertices_(std::move(colorVertices)),
          squareVertices_(std::move(squareVertices))
      {}

      void operator() (std::ostream& os) const {
        os << "  graph [fontname = \"Arial\", layout=\"dot\"];\n"
          << "  node [fontname = \"Arial\", shape = circle, style = filled];\n"
          << "  edge [fontname = \"Arial\"];\n"
          << R"(  labelloc="t"; label=")" << title_ << "\"\n";
      }

      void operator() (
        std::ostream& os,
        const typename DisplayGraph::vertex_descriptor& vertexIndex
      ) const {
        if(vertexIndex == rootIndex) {
          os << R"([label=")" << title_ << "\\n-\\n" << graphBase_[rootIndex].representation << R"(")";
        } else {
          os << R"([label=")" << graphBase_[vertexIndex].representation << R"(")";
        }

        if(colorVertices_.count(vertexIndex) > 0) {
          os << R"(, fillcolor="tomato", fontcolor="white")";
        }

        if(squareVertices_.count(vertexIndex) > 0) {
          os << R"(, shape="square")";
        }

        os << R"(])";
      }

      void operator() (
        std::ostream& /* os */,
        const typename DisplayGraph::edge_descriptor& /* edge */
      ) const { /* do nothing particular for edges */ }

    };

    SmallGraphWriter propertyWriter {graph, title, colorVertices, squareVertices};

    std::stringstream ss;

    boost::write_graphviz(
      ss,
      graph,
      propertyWriter,
      propertyWriter,
      propertyWriter
    );

    return ss.str();
  }

  std::string _make4BGraph(
    const TreeVertexIndex& sourceIndex,
    const std::map<
      TreeVertexIndex,
      std::set<VariantType>
    >& representativeStereodescriptors,
    const TreeVertexIndex& branchA,
    const TreeVertexIndex& branchB,
    const std::vector<
      std::vector<VariantType>
    >& branchAOrders,
    const std::vector<
      std::vector<VariantType>
    >& branchBOrders,
    const std::vector<
      std::vector<VariantType>
    >::reverse_iterator& branchAIter,
    const std::vector<
      std::vector<VariantType>
    >::reverse_iterator& branchBIter
  ) const;

  std::vector<
    std::vector<AtomIndex>
  > mapToAtomIndices_(
    const std::vector<
      std::vector<TreeVertexIndex>
    >& treeRankingSets
  ) const;

  static bool relevantSeeds_(
    const std::map<
      TreeVertexIndex,
      std::vector<TreeVertexIndex>
    >& seeds,
    const std::vector<
      std::vector<TreeVertexIndex>
    >& undecidedSets
  );

  void applySequenceRules_(
    const boost::optional<AngstromPositions>& positionsOption
  );

  std::vector<
    std::vector<TreeVertexIndex>
  > auxiliaryApplySequenceRules_(
    TreeVertexIndex sourceIndex,
    const std::vector<TreeVertexIndex>& adjacentsToRank,
    const boost::optional<unsigned>& depthLimitOptional = boost::none
  ) const;

public:

  RankingTree(
    const Graph& graph,
    const StereopermutatorList& stereopermutators,
    std::string molGraphviz,
    AtomIndex atomToRank,
    const std::vector<AtomIndex>& excludeIndices = {},
    ExpansionOption expansionMethod = ExpansionOption::OnlyRequiredBranches,
    const boost::optional<AngstromPositions>& positionsOption = boost::none
  );


  std::vector<
    std::vector<AtomIndex>
  > getRanked() const;

  std::string dumpGraphviz(
    const std::string& title = "",
    const std::unordered_set<TreeVertexIndex>& squareVertices = {},
    const std::unordered_set<TreeVertexIndex>& colorVertices = {},
    const std::set<TreeEdgeIndex>& colorEdges = {}
  ) const;
};

} // namespace Molassembler
} // namespace Scine

#endif