diff --git a/bnb/edmonds_karp.cpp b/bnb/edmonds_karp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1252cf5a2b6737cac834d027baaf98f51bb481fc
--- /dev/null
+++ b/bnb/edmonds_karp.cpp
@@ -0,0 +1,135 @@
+#include <limits>
+#include <queue>
+#include <cassert>
+
+#include <gp-bnb/edmonds_karp.hpp>
+
+edmonds_karp::edmonds_karp(const graph &g, node_id source, node_id sink)
+ : g(g), source(source), sink(sink) {
+ assert(source > 0); // Nodes have ids: 1, .., n
+ assert(source < g.num_nodes());
+};
+
+/* Indexes edges of the graph
+Every edge in the graph is mapped to its unique id
+indexed_edges std::map consists of key: pair<node_id, node_id> value: edge_id */
+void edmonds_karp::index_edges() {
+
+ unsigned int num_nodes = g.num_nodes();
+ node_id u = 1;
+ unsigned int eid = 0;
+
+ for (unsigned int i = 0; i < num_nodes; i++) {
+
+ std::vector<node_id> neighbors = g.get_adjacency(u);
+
+ for (node_id v : neighbors) {
+ auto node_pair = std::make_pair(u, v);
+
+ // edge was already listed
+ if (indexed_edges.count(std::make_pair(v, u)) == 1) {
+ indexed_edges[node_pair] = indexed_edges.at(std::make_pair(v, u));
+
+ } else {
+ indexed_edges[node_pair] = eid;
+ eid++;
+ }
+ }
+ u++;
+ }
+}
+
+/* Breadth-first search in the graph from source to sink
+When the sink is reached gain is added to gain[sink] */
+int edmonds_karp::bfs(std::vector<int> &resid_flow, std::vector<unsigned int> &pred) const {
+
+ pred.clear();
+ pred.resize(g.num_nodes() + 1, -1 ); // undiscovered nodes are marked with -1
+ std::vector<int> gain(g.num_nodes(), 0.0);
+
+ std::queue<node_id> q;
+ q.push(source);
+ pred[source] = source;
+ gain[source] = std::numeric_limits<int>::max();
+ while (!q.empty()) {
+ node_id u = q.front();
+ q.pop();
+
+ bool sink_reached = false;
+
+ std::vector<node_id> neighbors = g.get_adjacency(u);
+
+ // iterate through neighbors of u
+ for (node_id v : neighbors) {
+ unsigned int edge_id = indexed_edges.at(std::make_pair(u, v));
+ int edge_weight = 1; // unweighted graph
+
+ if (((u >= v && flow[edge_id] < edge_weight)
+ || ( u < v && resid_flow[edge_id] < edge_weight))
+ && pred[v] == (unsigned int) -1) { // only add those neighbors with rest capacity and which were not discovered yet
+ pred[v] = u;
+ gain[v] = std::min(gain[u], edge_weight - (u >= v ? flow[edge_id] : resid_flow[edge_id]));
+ if (v != sink && !sink_reached) {
+ q.push(v);
+ } else {
+ sink_reached = true;
+ }
+ }
+ }
+
+ if (sink_reached) {
+ return gain[sink];
+ }
+ }
+ return 0.0;
+};
+
+
+/* Edmonds-Karp Algorithm for unweighted, undirected graphs
+Step 0: Index graph edges
+Step 1: Initialize flow and residual flow vectors with length |E|
+--Loop:
+Step 2: Perform BFS that returns max gain from source to sink in one path
+Step 3: Add gain that was calculated during BFS to the flow value. Break from the loop if gain was 0
+Step 4: Update flow and residual flow values for each node */
+void edmonds_karp::run() {
+ index_edges();
+ int num_edges = indexed_edges.size() / 2;
+
+ flow.clear();
+ flow.resize(num_edges, 0.0);
+ std::vector<int> resid_flow(num_edges, 0.0);
+
+ flow_value = 0;
+
+ while (true) {
+ std::vector<node_id> pred;
+ // Perform BFS that returns max gain from source to sink in one path
+ int gain = bfs(resid_flow, pred);
+
+ // If gain == 0.0 that means there exist no more paths from source to sink
+ if (gain == 0) break;
+
+ // Add gain that was calculated during BFS to the flow value
+ flow_value += gain;
+ node_id v = sink;
+
+ // Update flow and residual flow values for each edge
+ while (v != source) {
+ node_id u = pred[v];
+ int edge_id = indexed_edges[std::make_pair(u, v)];
+ if (u >= v) {
+ flow[edge_id] += gain;
+ resid_flow[edge_id] -= gain;
+ } else {
+ flow[edge_id] -= gain;
+ resid_flow[edge_id] += gain;
+ }
+ v = u;
+ }
+ }
+};
+
+int edmonds_karp::get_max_flow() const {
+ return flow_value;
+};
\ No newline at end of file
diff --git a/include/gp-bnb/edmonds_karp.hpp b/include/gp-bnb/edmonds_karp.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..2c1fac480d026ab921af78bd3151c7ae176d243d
--- /dev/null
+++ b/include/gp-bnb/edmonds_karp.hpp
@@ -0,0 +1,56 @@
+#ifndef EDMONDSKARP_H_
+#define EDMONDSKARP_H_
+
+#include "graph.hpp"
+#include <vector>
+#include <map>
+
+class edmonds_karp {
+private:
+ const graph &g;
+
+ std::map<std::pair<node_id, node_id>, int> indexed_edges;
+
+ node_id source;
+ node_id sink;
+
+ int flow_value;
+ std::vector<int> flow;
+
+ /**
+ * Indexes edges of the graph
+ * Every edge in the graph is mapped to its unique id
+ */
+ void index_edges();
+ /**
+ * Performs a breadth-first search on the graph from the source node to find an augmenting path to the sink node respecting the flow values
+ * @param residFlow The residual flow in the network.
+ * @param pred Used to store the path from the source to the sink.
+ * @return The gain in terms of flow.
+ */
+ int bfs(std::vector<int> &resid_flow, std::vector<unsigned int> &pred) const;
+
+public:
+ /**
+ * Constructs an instance of the EdmondsKarp algorithm for the given graph, source and sink
+ * @param graph The graph.
+ * @param source The source node.
+ * @param sink The sink node.
+ */
+ edmonds_karp(const graph &g, node_id source, node_id sink);
+
+ /**
+ * Computes the maximum flow, executes the EdmondsKarp algorithm.
+ * For unweighted, undirected Graphs.
+ */
+ void run();
+
+ /**
+ * Returns the value of the maximum flow from source to sink.
+ *
+ * @return The maximum flow value
+ */
+ int get_max_flow() const;
+};
+
+#endif /* EDMONDSKARP_H_ */
\ No newline at end of file