summaryrefslogtreecommitdiff
path: root/src/analysis_tools.cpp
blob: dea20f0167eb735b46f8f6d9df8e76feba8c21a7 (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
124
125
126
127
128
129
130
131

#include "analysis_tools.hpp"

class badcycleException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "Could not find a valid cycle on the broken system.";
  }
} badcycleex;

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<std::pair<std::array<unsigned, 2>, 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);
    }
  }

  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);
    }
  }

  std::list<std::pair<std::array<unsigned, 2>, std::list<unsigned>>> output;

  // 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) {
      output.push_back({{1, 0}, cycle});
    } else if (crossing_count[0] % 2 == 0 && crossing_count[1] % 2 == 1) {
      output.push_back({{0, 1}, cycle});
    }
  }

  return output;
}