summaryrefslogtreecommitdiff
path: root/src/analysis_tools.cpp
blob: 778e713458a5312709230934fe7d004c96f9c51f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

#include "analysis_tools.hpp"

template <class T>
bool is_shorter(const std::list<T> &l1, const std::list<T> &l2) {
  return l1.size() < l2.size();
}

bool trivial(boost::detail::edge_desc_impl<boost::undirected_tag,unsigned long>) {
  return true;
}

std::list<unsigned> find_minimal_crack(const Graph& G, const network& n) {
  Graph Gtmp(n.G.vertices.size());
  std::list<unsigned> removed_edges;

  class add_tree_edges : public boost::default_dfs_visitor {
    public:
      Graph& G;
      std::list<unsigned>& E;

      add_tree_edges(Graph& G, std::list<unsigned>& E) : G(G), E(E) {}

      void tree_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        boost::add_edge(boost::source(e, g), boost::target(e, g), g[e], G);
      }

      void back_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        if (!(boost::edge(boost::source(e, g), boost::target(e, g), G).second)) {
          E.push_back(g[e].index);
        }
      }
  };

  add_tree_edges ate(Gtmp, removed_edges);
  boost::depth_first_search(G, visitor(ate));

  class find_cycle : public boost::default_dfs_visitor {
    public:
      std::list<unsigned>& E;
      unsigned end;
      struct done{};

      find_cycle(std::list<unsigned>& E, unsigned end) : E(E), end(end) {}

      void discover_vertex(boost::graph_traits<Graph>::vertex_descriptor v, const Graph& g) {
        if (v == end) {
          throw done{};
        }
      }

      void examine_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        E.push_back(g[e].index);
      }

      void finish_edge(boost::graph_traits<Graph>::edge_descriptor e, const Graph& g) {
        E.erase(std::find(E.begin(), E.end(), g[e].index));
      }
  };

  std::list<std::list<unsigned>> cycles;

  for (auto edge : removed_edges) {
    std::list<unsigned> cycle = {edge};
    find_cycle vis(cycle, n.G.dual_edges[edge].v[1]);
    std::vector<boost::default_color_type> new_color_map(boost::num_vertices(Gtmp));
    try {
    boost::depth_first_visit(Gtmp, n.G.dual_edges[edge].v[0], vis, boost::make_iterator_property_map(new_color_map.begin(), boost::get(boost::vertex_index, Gtmp), new_color_map[0]));
    } catch(find_cycle::done const&) {
      cycles.push_back(cycle);
    }
  }

  if (cycles.size() > 1) {
    std::list<std::valarray<uint8_t>> bool_cycles;
    for (auto cycle : cycles) {
      std::valarray<uint8_t> bool_cycle(n.G.edges.size());
      for (auto v : cycle) {
        bool_cycle[v] = 1;
      }
      bool_cycles.push_back(bool_cycle);
    }

    // generate all possible cycles by taking xor of the edge sets of the known cycles
    for (auto it1 = bool_cycles.begin(); it1 != std::prev(bool_cycles.end()); it1++) {
      for (auto it2 = std::next(it1); it2 != bool_cycles.end(); it2++) {
        std::valarray<uint8_t> new_bool_cycle = (*it1) ^ (*it2);
        std::list<unsigned> new_cycle;
        unsigned pos = 0;
        for (uint8_t included : new_bool_cycle) {
          if (included) {
            new_cycle.push_back(pos);
          }
          pos++;
        }
        cycles.push_back(new_cycle);
      }
    }

    // find the cycle representing the crack by counting boundary crossings
    for (auto cycle : cycles) {
      std::array<unsigned, 2> crossing_count{0,0};

      for (auto edge : cycle) {
        if (n.G.dual_edges[edge].crossings[0]) {
          crossing_count[0]++;
        }
        if (n.G.dual_edges[edge].crossings[1]) {
          crossing_count[1]++;
        }
      }

      if (crossing_count[0] % 2 == 1 && crossing_count[1] % 2 == 0) {
        return cycle;
      }
    }
  } else {
    return cycles.front();
  }

  exit(5);
}