diff --git a/agraphbench2/makefile b/agraphbench2/makefile
new file mode 100644
index 0000000000000000000000000000000000000000..f873a5730886517043259346b3a0318c5dd3d6f0
--- /dev/null
+++ b/agraphbench2/makefile
@@ -0,0 +1,160 @@
+SHELL:=/bin/bash
+USE_RAMDISK ?= 0
+-include makefile.conf
+
+define NEW_LINE
+
+
+endef
+
+export NEW_LINE
+
+LDFLAGS_DEBUG:=$(addprefix -L, $(addsuffix lib-debug, $(LINK_PATHS))) -rdynamic $(addprefix -l, $(LINK_LIBRARIES)) -Wl,-rpath,.
+
+LDFLAGS_RELEASE:=$(addprefix -L, $(addsuffix lib-release, $(LINK_PATHS))) -rdynamic $(addprefix -l, $(LINK_LIBRARIES)) -Wl,-rpath,.
+
+OBJECTS_DEBUG:=$(patsubst src/%.cpp, obj-debug/%.o, $(shell find src/ -name *cpp))
+
+OBJECTS_RELEASE:=$(patsubst src/%.cpp, obj-release/%.o, $(shell find src/ -name *cpp))
+
+.PHONY: all build-debug clean-debug doc
+
+all:
+	@echo "What to do master?"
+
+obj%/makefile: makefile makefile.conf
+	if [ ! -w $(dir $@) ] && [ $(USE_RAMDISK) -eq 1 ]; then\
+		ln -s /tmp/`date +'%s%N'`-$(dir $@) $(subst /, , $(dir $@)) 2>/dev/null;\
+	fi;\
+	if [ -L $(subst /, , $(dir $@)) ]; then\
+		mkdir -p `readlink $(subst /, , $(dir $@))`;\
+	else\
+		mkdir -p $(dir $@);\
+	fi
+	echo "\
+	SHELL:=/bin/bash$${NEW_LINE}\
+	SRCDIR:=$${NEW_LINE}\
+	$${NEW_LINE}\
+	define NEW_LINE$${NEW_LINE}\
+	$${NEW_LINE}\
+	$${NEW_LINE}\
+	endef$${NEW_LINE}\
+	$${NEW_LINE}\
+	export NEW_LINE$${NEW_LINE}\
+	$${NEW_LINE}\
+	CXXFLAGS:= -pipe -std=c++11 \$$(CXX_OTHER_FLAGS) -c -Wall -pedantic -Wextra -Werror -fPIC $(addprefix -I, $(INCLUDE_PATHS))$${NEW_LINE}\
+	$${NEW_LINE}\
+	SOURCES:= \$$(shell find \$$(SOURCES_BASE_DIR)/\$$(SRCDIR) -maxdepth 1 -type f -name \"*.cpp\")$${NEW_LINE}\
+	DEPENDENCIES:= \$$(patsubst \$$(SOURCES_BASE_DIR)/\$$(SRCDIR)%.cpp, \$$(OBJECTS_BASE_DIR)/\$$(SRCDIR)%.d, \$$(SOURCES))$${NEW_LINE}\
+	OBJECTS:= \$$(patsubst %.d, %.o, \$$(DEPENDENCIES))$${NEW_LINE}\
+	SOURCES_DIRS:= \$$(shell find \$$(SOURCES_BASE_DIR)/\$$(SRCDIR) -maxdepth 1 -mindepth 1 -type d)$${NEW_LINE}\
+	OBJECTS_DIRS:= \$$(patsubst \$$(SOURCES_BASE_DIR)/\$$(SRCDIR)%, %/, \$$(SOURCES_DIRS))$${NEW_LINE}\
+	OBJECTS_DIRS_MAKEFILES:= \$$(patsubst %, %makefile, \$$(OBJECTS_DIRS))$${NEW_LINE}\
+	$${NEW_LINE}\
+	.PHONY: all$${NEW_LINE}\
+	.PRECIOUS: \$$(DEPENDECIES) \$$(OBJECTS_DIRS_MAKEFILES)$${NEW_LINE}\
+	$${NEW_LINE}\
+	all: \$$(OBJECTS_DIRS) \$$(OBJECTS)$${NEW_LINE}\
+	$${NEW_LINE}\
+	%.d: makefile$${NEW_LINE}\
+		@echo \"\\$${NEW_LINE}\
+		\$$(shell sha1sum <<< \"\$$@\" | sed \"s/  -//g\") = \\$$\$$(shell (\\$$\$$(CXX) -M \\$$\$$(CXXFLAGS) \$$(patsubst \$$(OBJECTS_BASE_DIR)/\$$(SRCDIR)%.d,\$$(SOURCES_BASE_DIR)/\$$(SRCDIR)%.cpp, \$$@) 2>/dev/null || echo \\\"\$$(patsubst \$$(OBJECTS_BASE_DIR)/\$$(SRCDIR)%.d,\$$(SOURCES_BASE_DIR)/\$$(SRCDIR)%.cpp, \$$@) FORCE\\\") | sed \\\"s/.*://g;s/\\\\\\\\\\\\\\\\//g\\\")\$$\$${NEW_LINE}\\$${NEW_LINE}\
+		\$$(patsubst %.d,%.o, \$$@): \\$$\$$(\$$(shell sha1sum <<< \"\$$@\" | sed \"s/  -//g\")) makefile\$$\$${NEW_LINE}\\$${NEW_LINE}\
+			\\$$\$$(CXX) \\$$\$$(CXXFLAGS) \\$$\$$< -o \$$(patsubst %.d,%.o, \$$@)\$$\$${NEW_LINE}\\$${NEW_LINE}\
+		\" > \$$@$${NEW_LINE}\
+	$${NEW_LINE}\
+	%/makefile: makefile$${NEW_LINE}\
+		mkdir -p \$$(dir \$$@)$${NEW_LINE}\
+		cp makefile \$$@$${NEW_LINE}\
+	$${NEW_LINE}\
+	%/: FORCE | %/makefile$${NEW_LINE}\
+		@accesstime=\`stat -c %Y \$$@\` && \\$${NEW_LINE}\
+		\$$(MAKE) -C \$$@ SRCDIR=\$$(SRCDIR)\$$(notdir \$$(patsubst %/, %, \$$@))/ OBJECTS_BASE_DIR=\$$(OBJECTS_BASE_DIR) SOURCES_BASE_DIR=\$$(SOURCES_BASE_DIR) CXX_OTHER_FLAGS=\"\$$(CXX_OTHER_FLAGS)\" && \\$${NEW_LINE}\
+		accesstime2=\`stat -c %Y \$$@\` && \\$${NEW_LINE}\
+		if [ "\$$\$$accesstime" -ne "\$$\$$accesstime2" ]; then \\$${NEW_LINE}\
+			touch .; \\$${NEW_LINE}\
+		fi$${NEW_LINE}\
+	$${NEW_LINE}\
+	FORCE:$${NEW_LINE}\
+	$${NEW_LINE}\
+	-include \$$(DEPENDENCIES)" > $@
+
+debug: build-debug
+
+release: build-release
+
+clean: clean-debug clean-release
+	$(RM) -r doc
+
+
+
+bin-debug/$(EXECUTABLE): obj-debug/ $(OBJECTS_DEBUG)
+	if [ ! -w $(dir $@) ] && [ $(USE_RAMDISK) -eq 1 ]; then\
+		ln -s /tmp/`date +'%s%N'`-$(dir $@) $(subst /, , $(dir $@)) 2>/dev/null;\
+	fi;\
+	if [ -L $(subst /, , $(dir $@)) ]; then\
+		mkdir -p `readlink $(subst /, , $(dir $@))`;\
+	else\
+		mkdir -p $(dir $@);\
+	fi
+	$(CXX) $(OBJECTS_DEBUG) -o $@ $(LDFLAGS_DEBUG)
+
+bin-release/$(EXECUTABLE): obj-release/ $(OBJECTS_RELEASE)
+	if [ ! -w $(dir $@) ] && [ $(USE_RAMDISK) -eq 1 ]; then\
+		ln -s /tmp/`date +'%s%N'`-$(dir $@) $(subst /, , $(dir $@)) 2>/dev/null;\
+	fi;\
+	if [ -L $(subst /, , $(dir $@)) ]; then\
+		mkdir -p `readlink $(subst /, , $(dir $@))`;\
+	else\
+		mkdir -p $(dir $@);\
+	fi
+	$(CXX) $(OBJECTS_RELEASE) -o $@ $(LDFLAGS_RELEASE)
+
+
+
+obj-debug/: FORCE | obj-debug/makefile
+	$(MAKE) -C $@ OBJECTS_BASE_DIR=$(realpath obj-debug) SOURCES_BASE_DIR=$(realpath src) CXX_OTHER_FLAGS="-g -O0 -DDEBUG"
+
+obj-release/: FORCE | obj-release/makefile
+	$(MAKE) -C $@ OBJECTS_BASE_DIR=$(realpath obj-release) SOURCES_BASE_DIR=$(realpath src) CXX_OTHER_FLAGS="-O3 -DNDEBUG -DRELEASE"
+
+
+
+$(OBJECTS_DEBUG): obj-debug/
+
+$(OBJECTS_RELEASE): obj-release/
+
+
+build-debug: bin-debug/$(EXECUTABLE)
+
+build-release: bin-release/$(EXECUTABLE)
+
+
+
+clean-debug:
+	if [ -L obj-debug ]; then\
+		rm -r `readlink obj-debug`;\
+	fi
+	if [ -L bin-debug ]; then\
+		rm -r `readlink bin-debug`;\
+	fi
+	$(RM) -r *.o *.d bin-debug obj-debug
+
+clean-release:
+	if [ -L obj-release ]; then\
+		rm -r `readlink obj-release`;\
+	fi
+	if [ -L lib-release ]; then\
+		rm -r `readlink lib-release`;\
+	fi
+	$(RM) -r *.o *.d bin-release obj-release
+
+
+
+FORCE:
+
+
+
+doc:
+	doxygen
+
diff --git a/agraphbench2/makefile.conf b/agraphbench2/makefile.conf
new file mode 100644
index 0000000000000000000000000000000000000000..9b09eee8f57d74a38d64c15ded3263fbb5ca7415
--- /dev/null
+++ b/agraphbench2/makefile.conf
@@ -0,0 +1,4 @@
+EXECUTABLE:=agraphbench2
+LINK_PATHS=../alib2elgo/ ../alib2algo/ ../alib2data/ ../alib2std/
+LINK_LIBRARIES=alib2elgo alib2algo alib2data alib2std xml2
+INCLUDE_PATHS=\$$(SOURCES_BASE_DIR)/../../alib2elgo/src/ \$$(SOURCES_BASE_DIR)/../../alib2algo/src/ \$$(SOURCES_BASE_DIR)/../../alib2data/src/ \$$(SOURCES_BASE_DIR)/../../alib2std/src/ /usr/include/libxml2/
diff --git a/agraphbench2/src/agraphbench.cpp b/agraphbench2/src/agraphbench.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8efb0e1c05df933f1616c9b02dca1a0fe2c0dee0
--- /dev/null
+++ b/agraphbench2/src/agraphbench.cpp
@@ -0,0 +1,65 @@
+#include <tclap/CmdLine.h>
+#include <common/GlobalData.h>
+#include <string>
+#include <exception/AlibException.h>
+#include <factory/XmlDataFactory.hpp>
+#include <sax/SaxParseInterface.h>
+#include <sax/ParserException.h>
+#include <graph/isomorphism/Isomorphism.h>
+#include <graph/isomorphism/HopcroftWong.h>
+#include <graph/embedding/HopcroftTarjan.h>
+#include <graph/generate/RandomGraphFactory.h>
+
+int main(int argc, char** argv) {
+	try {
+		TCLAP::CmdLine cmd("Graph benchmark binary", ' ', "0.01");
+
+		TCLAP::ValueArg<int> size(	"",	"size",	"Size of graph",		false,	10, "integer");
+		cmd.add( size );
+
+		TCLAP::SwitchArg onlyhopcroft(		"",	"onlyhopcroft",	"Use only HopcroftWong algorithm",			false);
+		cmd.add( onlyhopcroft );
+
+		TCLAP::SwitchArg verbose(		"v",	"verbose",	"Be verbose",			false);
+		cmd.add( verbose );
+
+		cmd.parse(argc, argv);
+
+		common::GlobalData::measure = true;
+		if(verbose.isSet())
+			common::GlobalData::verbose = true;
+
+		graph::UndirectedGraph ug1 = graph::generate::RandomGraphFactory::generateUndirectedTriConnectedPlanarGraph( size.getValue() );
+		ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+		graph::UndirectedGraph ug2 = graph::generate::RandomGraphFactory::generateUndirectedIsomorphicGraph(ug1);
+		ug2.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug2));
+
+		std::chrono::measurements::start("HopcroftWong", std::chrono::measurements::Type::OVERALL);
+		if (!graph::isomorphism::HopcroftWong::hopcroftwong(ug1, ug2))
+			std::cout << "FAIL: hopcroftwong returned false!" << std::endl;
+		std::chrono::measurements::end();
+		std::clog << std::chrono::measurements::results() << std::endl;
+
+		if (!onlyhopcroft.isSet()) {
+			std::chrono::measurements::start("RecursiveByDefinition", std::chrono::measurements::Type::OVERALL);
+			if (!graph::isomorphism::Isomorphism::isomorphism(ug1, ug2))
+				std::cout << "FAIL: isomorphism returned false!" << std::endl;
+			std::chrono::measurements::end();
+			std::clog << std::chrono::measurements::results() << std::endl;
+		}
+
+		return 0;
+	} catch(const exception::AlibException& exception) {
+		alib::XmlDataFactory::toStdout(exception);
+		return 1;
+	} catch(const TCLAP::ArgException& exception) {
+		std::cerr << exception.error() << std::endl;
+		return 2;
+	} catch (const std::exception& exception) {
+		std::cerr << "Exception caught: " << exception.what() << std::endl;
+		return 3;
+	} catch(...) {
+		std::cerr << "Unknown exception caught." << std::endl;
+		return 127;
+	}
+}
diff --git a/alib2algo/src/graph/embedding/HopcroftTarjan.cpp b/alib2algo/src/graph/embedding/HopcroftTarjan.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ab61a0d0f1f032cc22e967c2d4c6aad1ae457fe4
--- /dev/null
+++ b/alib2algo/src/graph/embedding/HopcroftTarjan.cpp
@@ -0,0 +1,675 @@
+#include "HopcroftTarjan.h"
+
+#include <list>
+#include <stack>
+#include <algorithm>
+
+#include <exception/AlibException.h>
+
+namespace graph
+{
+
+namespace embedding
+{
+
+// Helpers
+
+using node = Node;
+using edge = DirectedEdge;
+
+template <typename T>
+using edge_array = std::unordered_map<edge, T>;
+
+template <typename T>
+using node_array = std::unordered_map<node, T>;
+
+template<typename T>
+inline void conc(T &container, T &l)
+{
+	container.splice(container.end(), l);
+}
+
+class GraphWrapper {
+public:
+	explicit GraphWrapper(DirectedGraph *g) : m_g(g) { }
+
+	std::vector<DirectedEdge> adj_edges(const node &n) const
+	{
+		std::vector<DirectedEdge> out;
+		const auto &edges = m_g->neighborEdges(n);
+		std::copy(edges.begin(), edges.end(), std::back_inserter(out));
+		sort(out);
+		return out;
+	}
+
+	std::vector<DirectedEdge> all_edges() const
+	{
+		std::vector<DirectedEdge> out;
+		const auto &edges = m_g->getEdges();
+		std::copy(edges.begin(), edges.end(), std::back_inserter(out));
+		sort(out);
+		return out;
+	}
+
+	edge first_adj_edge(const node &n) const
+	{
+		return adj_edges(n).front();
+	}
+
+	std::set<Node> all_nodes() const
+	{
+		return m_g->getNodes();
+	}
+
+	node first_node() const
+	{
+		return *m_g->getNodes().begin();
+	}
+
+	int number_of_nodes() const
+	{
+		return m_g->getNodes().size();
+	}
+
+	int number_of_edges() const
+	{
+		return m_g->getEdges().size();
+	}
+
+	void sort_edges(const edge_array<int> &sorted)
+	{
+		m_sortedEdges = sorted;
+	}
+
+	void del_edge(const edge &e)
+	{
+		m_g->removeEdge(e);
+	}
+
+private:
+	// FIXME: This kills performance!
+	// The correct way to implement this is to use adjacency lists
+	// in graph wrapper and not a DirectedGraph.
+	void sort(std::vector<DirectedEdge> &edges) const
+	{
+		if (m_sortedEdges.empty())
+			return;
+
+		std::sort(edges.begin(), edges.end(), [this](const DirectedEdge &a, const DirectedEdge &b) {
+			return m_sortedEdges.at(a) < m_sortedEdges.at(b);
+		});
+	}
+
+	DirectedGraph *m_g;
+	edge_array<int> m_sortedEdges;
+};
+
+using graph = GraphWrapper;
+
+bool createReversal(const graph &g, edge_array<edge> &reversal)
+{
+	std::vector<DirectedEdge> edges = g.all_edges();
+
+	while (!edges.empty()) {
+		DirectedEdge edge = edges.back();
+		edges.pop_back();
+		bool found = false;
+		for (size_t i = 0; i < edges.size(); ++i) {
+			const DirectedEdge &e = edges[i];
+			if (edge.getFromNode() == e.getToNode() && edge.getToNode() == e.getFromNode()) {
+				reversal.insert({edge, e});
+				reversal.insert({e, edge});
+				edges.erase(edges.begin() + i);
+				found = true;
+				break;
+			}
+		}
+		if (!found)
+			return false;
+	}
+
+	return true;
+}
+
+//  Mehlhorn, Kurt and Mutzel, Petra and Näher, Stefan (1994)
+//  An Implementation of the Hopcroft and Tarjan Planarity Test and Embedding Algorithm.
+
+const int LEFT = 1;
+const int RIGHT = 2;
+
+class Block
+{
+private:
+	std::list<int> Latt, Ratt; // list of attachments
+	std::list<edge> Lseg, Rseg; // list of segments represented by their defining edges
+
+public:
+	Block(const edge &e, std::list<int> &A)
+	{
+		Lseg.push_back(e);
+		conc(Latt, A);
+		// the other two lists are empty
+	}
+
+	void flip()
+	{
+		std::list<int> ha;
+		std::list<edge> he;
+		// we first interchange |Latt| and |Ratt| and then |Lseg| and |Rseg|
+		conc(ha, Ratt);
+		conc(Ratt, Latt);
+		conc(Latt, ha);
+		conc(he, Rseg);
+		conc(Rseg, Lseg);
+		conc(Lseg, he);
+	}
+
+	int head_of_Latt()
+	{
+		return Latt.front();
+	}
+
+	bool empty_Latt()
+	{
+		return Latt.empty();
+	}
+
+	int head_of_Ratt()
+	{
+		return Ratt.front();
+	}
+
+	bool empty_Ratt()
+	{
+		return Ratt.empty();
+	}
+
+	bool left_interlace(std::stack<Block*> &S)
+	{
+		// check for interlacing with the left side of the topmost block of |S|
+		if (Latt.empty())
+			throw exception::AlibException("HopcroftTarjan: Latt is never empty");
+
+		if (!S.empty() && !((S.top())->empty_Latt()) && Latt.back() < (S.top())->head_of_Latt())
+			return true;
+		else
+			return false;
+	}
+
+	bool right_interlace(std::stack<Block*> &S)
+	{
+		// check for interlacing with the right side of the topmost block of |S|
+		if (Latt.empty())
+			throw exception::AlibException("HopcroftTarjan: Latt is never empty");
+
+		if (!S.empty() && !((S.top())->empty_Ratt()) && Latt.back() < (S.top())->head_of_Ratt())
+			return true;
+		else
+			return false;
+	}
+
+	void combine(Block *Bprime)
+	{
+		// add block Bprime to the rear of |this| block
+		conc(Latt, Bprime->Latt);
+		conc(Ratt, Bprime->Ratt);
+		conc(Lseg, Bprime->Lseg);
+		conc(Rseg, Bprime->Rseg);
+		delete Bprime;
+	}
+
+	bool clean(int dfsnum_w, edge_array<int> &alpha)
+	{
+		// remove all attachments to |w|; there may be several
+		while (!Latt.empty() && Latt.front() == dfsnum_w)
+			Latt.pop_front();
+		while (!Ratt.empty() && Ratt.front() == dfsnum_w)
+			Ratt.pop_front();
+
+		if (!Latt.empty() || !Ratt.empty())
+			return false;
+
+		// |Latt| and |Ratt| are empty; we record the placement of the subsegments in |alpha|.
+		for (const edge &e : Lseg)
+			alpha[e] = LEFT;
+
+		for (const edge &e : Rseg)
+			alpha[e] = RIGHT;
+
+		return true;
+	}
+
+	void add_to_Att(std::list<int> &Att, int dfsnum_w0, edge_array<int> &alpha)
+	{
+		// add the block to the rear of |Att|. Flip if necessary
+		if (!Ratt.empty() && head_of_Ratt() > dfsnum_w0)
+			flip();
+
+		conc(Att, Latt);
+		conc(Att, Ratt);
+
+		// This needs some explanation. Note that |Ratt| is either empty
+		// or $\{w0\}$. Also if |Ratt| is non-empty then all subsequent sets are contained
+		// in $\{w0\}$. So we indeed compute an ordered set of attachments.
+
+		for (const edge &e : Lseg)
+			alpha[e] = LEFT;
+
+		for (const edge &e : Rseg)
+			alpha[e] = RIGHT;
+	}
+};
+
+static void dfs_in_reorder(const graph &G, std::list<edge> &Del, node v, int &dfs_count,
+	node_array<bool> &reached,
+	node_array<int> &dfsnum,
+	node_array<int> &lowpt1, node_array<int> &lowpt2,
+	node_array<node> &parent);
+
+static void reorder(graph &G, node_array<int> &dfsnum, node_array<node> &parent)
+{
+	node_array<bool> reached;
+	int dfs_count = 0;
+	std::list<edge> Del;
+	node_array<int> lowpt1, lowpt2;
+
+	dfs_in_reorder(G, Del, G.first_node(), dfs_count, reached, dfsnum, lowpt1, lowpt2, parent);
+
+	// remove forward and reversals of tree edges
+
+	for (const edge &e : Del)
+		G.del_edge(e);
+
+	// we now reorder adjacency lists as described in \cite[page 101]{Me:book}
+
+	edge_array<int> cost;
+	for (const edge &e : G.all_edges()) {
+		node v = e.getFromNode();
+		node w = e.getToNode();
+		cost[e] = ((dfsnum[w] < dfsnum[v]) ?
+		2 * dfsnum[w] :
+		((lowpt2[w] >= dfsnum[v]) ?
+		2 * lowpt1[w] : 2 * lowpt1[w] + 1));
+	}
+
+	G.sort_edges(cost);
+}
+
+static void dfs_in_reorder(const graph &G, std::list<edge> &Del, node v, int &dfs_count,
+	node_array<bool> &reached,
+	node_array<int> &dfsnum,
+	node_array<int> &lowpt1, node_array<int> &lowpt2,
+	node_array<node> &parent)
+{
+	dfsnum[v] = dfs_count++;
+	lowpt1[v] = lowpt2[v] = dfsnum[v];
+	reached[v] = true;
+
+	for (const edge &e : G.adj_edges(v)) {
+		node w = e.getToNode();
+		if (!reached[w]) { // e is a tree edge
+			parent.erase(w);
+			parent.insert({w, v});
+			dfs_in_reorder(G, Del, w, dfs_count, reached, dfsnum, lowpt1, lowpt2, parent);
+			lowpt1[v] = std::min(lowpt1[v], lowpt1[w]);
+		} // end tree edge
+		else {
+			lowpt1[v] = std::min(lowpt1[v], dfsnum[w]); // no effect for forward edges
+			if ((dfsnum[w] >= dfsnum[v]) || w == parent.at(v)) // forward edge or reversal of tree edge
+				Del.push_back(e);
+		} // end non-tree edge
+	} // end forall
+
+	// we know |lowpt1[v]| at this point and now make a second pass over all
+	// adjacent edges of |v| to compute |lowpt2|
+
+	for (const edge &e : G.adj_edges(v)) {
+		node w = e.getToNode();
+		if (parent.at(w) == v) { // tree edge
+			if (lowpt1[w] != lowpt1[v])
+				lowpt2[v] = std::min(lowpt2[v], lowpt1[w]);
+			lowpt2[v] = std::min(lowpt2[v], lowpt2[w]);
+		} // end tree edge
+		else if (lowpt1[v] != dfsnum[w]) // all other edges
+			lowpt2[v] = std::min(lowpt2[v], dfsnum[w]);
+	}
+}
+
+static bool strongly_planar(edge e0, const graph &G, std::list<int> &Att,
+	edge_array<int> &alpha,
+	node_array<int> &dfsnum,
+	node_array<node> &parent)
+{
+	node x = e0.getFromNode();
+
+	node y = e0.getToNode();
+
+	edge e = G.first_adj_edge(y);
+
+	node wk = y;
+
+	while (dfsnum[e.getToNode()] > dfsnum[wk]) { // e is a tree edge
+		wk = e.getToNode();
+		e = G.first_adj_edge(wk);
+	}
+
+	node w0 = e.getToNode();
+
+	node w = wk;
+
+	std::stack<Block*> S;
+
+	while (w != x) {
+		int count = 0;
+		for (const edge &e : G.adj_edges(w)) {
+			count++;
+			if (count != 1) { // no action for first edge
+				std::list<int> A;
+
+				if (dfsnum[w] < dfsnum[e.getToNode()]) { // tree edge
+					if (!strongly_planar(e, G, A, alpha, dfsnum, parent)) {
+						while (!S.empty()) {
+							delete S.top();
+							S.pop();
+						}
+						return false;
+					}
+				}
+				else {
+					A.push_back(dfsnum[e.getToNode()]); // a back edge
+				}
+
+				Block *B = new Block(e, A);
+
+				while (true) {
+					if (B->left_interlace(S))
+						(S.top())->flip();
+
+					if (B->left_interlace(S)) {
+						delete B;
+						while (!S.empty()) {
+							delete S.top();
+							S.pop();
+						}
+						return false;
+					};
+
+					if (B->right_interlace(S)) {
+						B->combine(S.top());
+						S.pop();
+					}
+					else
+						break;
+				} // end while
+
+				S.push(B);
+			} // end if
+		} // end forall
+
+		while (!S.empty() && (S.top())->clean(dfsnum[parent.at(w)], alpha)) {
+			delete S.top();
+			S.pop();
+		}
+
+		w = parent.at(w);
+	} // end while
+
+	Att.clear();
+	while (!S.empty()) {
+		Block *B = S.top();
+		S.pop();
+
+		if (!(B->empty_Latt()) && !(B->empty_Ratt()) &&
+		(B->head_of_Latt() > dfsnum[w0]) && (B->head_of_Ratt() > dfsnum[w0])) {
+			delete B;
+			while (!S.empty()) {
+				delete S.top();
+				S.pop();
+			}
+			return false;
+		}
+
+		B->add_to_Att(Att, dfsnum[w0], alpha);
+		delete B;
+	} // end while
+
+	// Let's not forget (as the book does) that $w0$ is an attachment of $S(e0)$
+	// except if $w0 = x$.
+
+	if (w0 != x)
+		Att.push_back(dfsnum[w0]);
+
+	return true;
+}
+
+static void embedding(edge e0, int t, const graph &G,
+	edge_array<int> &alpha,
+	node_array<int> &dfsnum,
+	std::list<edge> &T, std::list<edge> &A, int &cur_nr,
+	edge_array<int> &sort_num, node_array<edge> &tree_edge_into,
+	node_array<node> &parent, edge_array<edge> &reversal)
+{
+	node x = e0.getFromNode();
+
+	node y = e0.getToNode();
+
+	tree_edge_into.erase(y);
+	tree_edge_into.insert({y, e0});
+
+	edge e = G.first_adj_edge(y);
+
+	node wk = y;
+
+	while (dfsnum[e.getToNode()] > dfsnum[wk]) { // e is a tree edge
+		wk = e.getToNode();
+		tree_edge_into.erase(wk);
+		tree_edge_into.insert({wk, e});
+		e = G.first_adj_edge(wk);
+	}
+
+	node w0 = e.getToNode();
+	edge back_edge_into_w0 = e;
+
+	node w = wk;
+
+	std::list<edge> Al, Ar;
+	std::list<edge> Tprime, Aprime;
+
+	T.clear();
+	T.push_back(e); // |e = (wk,w0)| at this point, line (2)
+
+	while (w != x) {
+		int count = 0;
+		for (const edge &e : G.adj_edges(w)) {
+			count++;
+			if (count != 1) { // no action for first edge
+				if (dfsnum[w] < dfsnum[e.getToNode()]) { // tree edge
+					int tprime = ((t == alpha[e]) ? LEFT : RIGHT);
+					embedding(e, tprime, G, alpha, dfsnum, Tprime, Aprime, cur_nr, sort_num, tree_edge_into, parent, reversal);
+				}
+				else { // back edge
+					Tprime.push_back(e); // $e$
+					Aprime.push_back(reversal.at(e)); // reversal of $e$
+				}
+
+				if (t == alpha[e]) {
+					conc(Tprime, T);
+					conc(T, Tprime); // $T = Tprime\ conc\ T$
+					conc(Al, Aprime); // $Al = Al\ conc\ Aprime$
+				}
+				else {
+					conc(T, Tprime); // $ T\ = T\ conc\ Tprime $
+					conc(Aprime, Ar);
+					conc(Ar, Aprime); // $ Ar\ = Aprime\ conc\ Ar$
+				}
+			} // end if
+		} // end forall
+
+		T.push_back(reversal.at(tree_edge_into.at(w))); // $(w_{j-1},w_j)$
+
+		for (const edge &e : T)
+			sort_num[e] = cur_nr++;
+
+		// |w|'s adjacency list is now computed; we clear |T| and prepare for the
+		// next iteration by moving all darts incident to |parent[w]| from |Al| and
+		// |Ar| to |T|. These darts are at the rear end of |Al| and at the front end
+		// of |Ar|
+
+		T.clear();
+
+		// |parent[w]| is in |G|, |Al.tail| in |H|
+		while (!Al.empty() && Al.back().getFromNode() == parent.at(w)) {
+			T.push_front(Al.back()); // Pop removes from the rear
+			Al.pop_back();
+		}
+
+		T.push_back(tree_edge_into.at(w)); // push would be equivalent
+
+		while (!Ar.empty() && Ar.front().getFromNode() == parent.at(w)) {
+			T.push_back(Ar.front()); // pop removes from the front
+			Ar.pop_front();
+		}
+
+		w = parent.at(w);
+	} // end while
+
+	A.clear();
+
+	conc(A, Ar);
+	A.push_back(reversal.at(back_edge_into_w0));
+	conc(A, Al);
+}
+
+// |Gin| is a directed graph. |planar| decides whether |Gin| is planar.
+// If it is and |embed == true| then it also computes a
+// combinatorial embedding of |Gin| by suitably reordering
+// its adjacency lists.
+// |Gin| must be bidirected in that case.
+bool PLANAR(graph &G, bool embed, HopcroftTarjan::Result &res)
+{
+	int n = G.number_of_nodes();
+
+	if (n <= 3)
+		return true;
+
+	if (G.number_of_edges() > 6 * n - 12)
+		return false;
+
+	{
+		node_array<int> nr;
+		edge_array<int> cost;
+		int cur_nr = 0;
+		int n = G.number_of_nodes();
+
+		for (const node &v : G.all_nodes())
+			nr[v] = cur_nr++;
+
+		for (const edge &e : G.all_edges())
+			cost[e] = ((nr[e.getFromNode()] < nr[e.getToNode()]) ?
+			n * nr[e.getFromNode()] + nr[e.getToNode()] :
+			n * nr[e.getToNode()] + nr[e.getFromNode()]);
+
+		G.sort_edges(cost);
+
+		std::vector<edge> L = G.all_edges();
+		while (!L.empty()) {
+			edge e = L.front();
+			L.erase(L.begin());
+			// check whether the first edge on |L| is equal to the reversal of |e|. If so,
+			// delete it from |L|, if not, add the reversal to |G|
+			if (!L.empty() && (e.getFromNode() == L.front().getToNode()) && (L.front().getFromNode() == e.getToNode()))
+				L.erase(L.begin());
+			else
+				throw exception::AlibException("HopcroftTarjan: Graph not bidirected");
+		}
+	}
+
+	// An undirected planar graph has at most $3n-6$ edges; a directed graph may have twice as many
+
+	edge_array<edge> reversal;
+	if (!createReversal(G, reversal))
+		throw exception::AlibException("HopcroftTarjan: Graph not bidirected");
+
+	node_array<int> dfsnum;
+	node_array<node> parent;
+	for (const node &n : G.all_nodes())
+		parent.insert({n, node("HopcroftTarjan_invalid")}); // no default constructor -,-
+
+	reorder(G, dfsnum, parent);
+
+	edge_array<int> alpha;
+	for (const edge &e : G.all_edges())
+		alpha[e] = 0;
+
+	{
+		std::list<int> Att;
+
+		alpha[G.first_adj_edge(G.first_node())] = LEFT;
+
+		if (!strongly_planar(G.first_adj_edge(G.first_node()), G, Att, alpha, dfsnum, parent))
+			return false;
+	}
+
+	if (embed) {
+		std::list<edge> T, A; // lists of edges of |H|
+
+		int cur_nr = 0;
+		edge_array<int> sort_num;
+		node_array<edge> tree_edge_into;
+
+		embedding(G.first_adj_edge(G.first_node()), LEFT, G, alpha, dfsnum, T, A,
+		cur_nr, sort_num, tree_edge_into, parent, reversal);
+
+		// |T| contains all edges incident to the first node except the cycle edge into it.
+		// That edge comprises |A|
+
+		conc(T, A);
+
+		for (const edge &e : T)
+			sort_num[e] = cur_nr++;
+
+		std::vector<const DirectedEdge*> edgs(sort_num.size());
+		for (const auto &p : sort_num)
+			edgs[p.second] = &p.first;
+
+		for (const DirectedEdge *e : edgs)
+			res[e->getFromNode()].push_back(e->getToNode());
+	}
+
+	return true;
+}
+
+HopcroftTarjan::Result HopcroftTarjan::hopcrofttarjan(const Graph &graph)
+{
+	return getInstance().dispatch(graph.getData());
+}
+
+HopcroftTarjan::Result HopcroftTarjan::hopcrofttarjan(const DirectedGraph &graph)
+{
+	DirectedGraph copy = graph;
+	GraphWrapper graphw(&copy);
+	Result res;
+	if (!PLANAR(graphw, true, res))
+		throw exception::AlibException("HopcroftTarjan: Graph not planar");
+	return res;
+}
+
+HopcroftTarjan::Result HopcroftTarjan::hopcrofttarjan(const UndirectedGraph &graph)
+{
+	DirectedGraph copy;
+	for (const UndirectedEdge &e : graph.getEdges()) {
+		DirectedEdge e1(e.getFirstNode(), e.getSecondNode(), e.getName());
+		DirectedEdge e2(e.getSecondNode(), e.getFirstNode(), e.getName());
+		copy.addEdge(e1);
+		copy.addEdge(e2);
+	}
+	GraphWrapper graphw(&copy);
+	Result res;
+	if (!PLANAR(graphw, true, res))
+		throw exception::AlibException("HopcroftTarjan: Graph not planar");
+	return res;
+}
+
+} // namespace embedding
+
+} // namespace graph
diff --git a/alib2algo/src/graph/embedding/HopcroftTarjan.h b/alib2algo/src/graph/embedding/HopcroftTarjan.h
new file mode 100644
index 0000000000000000000000000000000000000000..797fdd6eed0efad3c731ba88a347c20e6b6cc8a7
--- /dev/null
+++ b/alib2algo/src/graph/embedding/HopcroftTarjan.h
@@ -0,0 +1,37 @@
+#ifndef GRAPH_HOPCROFT_TARJAN_H_
+#define GRAPH_HOPCROFT_TARJAN_H_
+
+#include <common/multipleDispatch.hpp>
+
+#include <graph/Graph.h>
+#include <graph/directed/DirectedGraph.h>
+#include <graph/undirected/UndirectedGraph.h>
+
+namespace graph
+{
+
+namespace embedding
+{
+
+// Computes combinatorial embedding of bidirected biconnected simple planar graph
+class HopcroftTarjan : public std::SingleDispatch<std::unordered_map<Node, std::vector<Node>>, graph::GraphBase>
+{
+public:
+	typedef std::unordered_map<Node, std::vector<Node>> Result;
+
+	static Result hopcrofttarjan(const Graph &graph);
+
+	static Result hopcrofttarjan(const DirectedGraph &graph);
+	static Result hopcrofttarjan(const UndirectedGraph &graph);
+
+	static HopcroftTarjan &getInstance() {
+		static HopcroftTarjan res;
+		return res;
+	}
+};
+
+} // namespace embedding
+
+} // namespace graph
+
+#endif // GRAPH_HOPCROFT_TARJAN_H_
diff --git a/alib2algo/src/graph/generate/RandomGraphFactory.cpp b/alib2algo/src/graph/generate/RandomGraphFactory.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0bb44e44b2173a03c7778e2f48d4af7bc54e4ed5
--- /dev/null
+++ b/alib2algo/src/graph/generate/RandomGraphFactory.cpp
@@ -0,0 +1,398 @@
+#include "RandomGraphFactory.h"
+
+#include <random>
+#include <algorithm>
+#include <unordered_map>
+#include <iostream>
+
+#ifndef NDEBUG
+#define GEN_GRAPH(x) std::cout << x << std::endl;
+#else
+#define GEN_GRAPH(x)
+#endif
+
+namespace graph
+{
+
+namespace generate
+{
+
+struct FuncOut
+{
+	UndirectedGraph graph;
+	std::vector<Node> targetNodes; // 3 nodes
+	std::vector<Node> sourceNodes; // 3 nodes
+};
+
+static int NODE_LABEL;
+
+// G1
+static FuncOut gen5_1()
+{
+	GEN_GRAPH("gen5_1");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n2, n3));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n5, n4));
+	g.addEdge(UndirectedEdge(n5, n3));
+	g.addEdge(UndirectedEdge(n5, n2));
+	return {g, {n2, n5, n4}, {n2, n1, n4}};
+}
+
+static FuncOut gen6_1()
+{
+	GEN_GRAPH("gen6_1");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n2, n3));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n5, n4));
+	g.addEdge(UndirectedEdge(n5, n3));
+	g.addEdge(UndirectedEdge(n5, n2));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n6, n1));
+	return {g, {n2, n5, n6}, {n2, n1, n6}};
+}
+
+static FuncOut gen7_1()
+{
+	GEN_GRAPH("gen7_1");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+	Node n7(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n2, n3));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n5, n4));
+	g.addEdge(UndirectedEdge(n5, n3));
+	g.addEdge(UndirectedEdge(n5, n2));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n6, n1));
+	g.addEdge(UndirectedEdge(n7, n2));
+	g.addEdge(UndirectedEdge(n7, n1));
+	g.addEdge(UndirectedEdge(n7, n6));
+	return {g, {n2, n5, n6}, {n2, n7, n6}};
+}
+
+static FuncOut gen8_1()
+{
+	GEN_GRAPH("gen8_1");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+	Node n7(NODE_LABEL++);
+	Node n8(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n2, n3));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n5, n4));
+	g.addEdge(UndirectedEdge(n5, n3));
+	g.addEdge(UndirectedEdge(n5, n2));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n6, n1));
+	g.addEdge(UndirectedEdge(n7, n2));
+	g.addEdge(UndirectedEdge(n7, n1));
+	g.addEdge(UndirectedEdge(n7, n6));
+	g.addEdge(UndirectedEdge(n8, n5));
+	g.addEdge(UndirectedEdge(n8, n2));
+	g.addEdge(UndirectedEdge(n8, n7));
+	return {g, {n8, n5, n6}, {n8, n7, n6}};
+}
+
+// G2
+static FuncOut gen5_2()
+{
+	GEN_GRAPH("gen5_2");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n5));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n4));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n5));
+	return {g, {n5, n1, n2}, {n5, n3, n4}};
+}
+
+static FuncOut gen6_2()
+{
+	GEN_GRAPH("gen6_2");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n5));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n4));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n3));
+	g.addEdge(UndirectedEdge(n6, n4));
+	return {g, {n5, n1, n2}, {n5, n6, n4}};
+}
+
+static FuncOut gen7_2()
+{
+	GEN_GRAPH("gen7_2");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+	Node n7(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n5));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n4));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n3));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n7, n2));
+	g.addEdge(UndirectedEdge(n7, n4));
+	g.addEdge(UndirectedEdge(n7, n6));
+	return {g, {n5, n1, n2}, {n5, n6, n7}};
+}
+
+static FuncOut gen8_2()
+{
+	GEN_GRAPH("gen8_2");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+	Node n7(NODE_LABEL++);
+	Node n8(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n1, n5));
+	g.addEdge(UndirectedEdge(n1, n3));
+	g.addEdge(UndirectedEdge(n1, n4));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n4));
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n6, n3));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n7, n2));
+	g.addEdge(UndirectedEdge(n7, n4));
+	g.addEdge(UndirectedEdge(n7, n6));
+	g.addEdge(UndirectedEdge(n8, n5));
+	g.addEdge(UndirectedEdge(n8, n1));
+	g.addEdge(UndirectedEdge(n8, n7));
+	return {g, {n5, n8, n7}, {n5, n6, n7}};
+}
+
+// G3
+static FuncOut gen5_3()
+{
+	GEN_GRAPH("gen5_3");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n1));
+	g.addEdge(UndirectedEdge(n3, n2));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n4, n1));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n5));
+	return {g, {n1, n4, n3}, {n1, n2, n5}};
+}
+
+static FuncOut gen6_3()
+{
+	GEN_GRAPH("gen6_3");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n1));
+	g.addEdge(UndirectedEdge(n3, n2));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n4, n1));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n5));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n6, n3));
+	g.addEdge(UndirectedEdge(n6, n5));
+	return {g, {n1, n4, n6}, {n2, n5, n6}};
+}
+
+static FuncOut gen7_3()
+{
+	GEN_GRAPH("gen7_3");
+
+	Node n1(NODE_LABEL++);
+	Node n2(NODE_LABEL++);
+	Node n3(NODE_LABEL++);
+	Node n4(NODE_LABEL++);
+	Node n5(NODE_LABEL++);
+	Node n6(NODE_LABEL++);
+	Node n7(NODE_LABEL++);
+
+	UndirectedGraph g;
+	g.addEdge(UndirectedEdge(n3, n4));
+	g.addEdge(UndirectedEdge(n3, n1));
+	g.addEdge(UndirectedEdge(n3, n2));
+	g.addEdge(UndirectedEdge(n3, n5));
+	g.addEdge(UndirectedEdge(n4, n1));
+	g.addEdge(UndirectedEdge(n1, n2));
+	g.addEdge(UndirectedEdge(n2, n5));
+	g.addEdge(UndirectedEdge(n6, n4));
+	g.addEdge(UndirectedEdge(n6, n3));
+	g.addEdge(UndirectedEdge(n6, n5));
+	g.addEdge(UndirectedEdge(n7, n1));
+	g.addEdge(UndirectedEdge(n7, n2));
+	g.addEdge(UndirectedEdge(n7, n6));
+	return {g, {n1, n4, n6}, {n1, n7, n6}};
+}
+
+static void mergeGraphs(FuncOut &target, const FuncOut &source)
+{
+	if (target.targetNodes.empty()) {
+		target = source;
+		return;
+	}
+
+	for (const UndirectedEdge &e : source.graph.getEdges())
+		target.graph.addEdge(e);
+
+	target.graph.addEdge(UndirectedEdge(target.targetNodes.at(0), source.sourceNodes.at(0)));
+	target.graph.addEdge(UndirectedEdge(target.targetNodes.at(1), source.sourceNodes.at(1)));
+	target.graph.addEdge(UndirectedEdge(target.targetNodes.at(2), source.sourceNodes.at(2)));
+	target.targetNodes = source.targetNodes;
+}
+
+UndirectedGraph RandomGraphFactory::generateUndirectedTriConnectedPlanarGraph(int size)
+{
+	NODE_LABEL = 0;
+
+	std::unordered_map<size_t, std::vector<std::function<FuncOut()>>> funcs;
+	funcs[5] = {gen5_1, gen5_2, gen5_3};
+	funcs[6] = {gen6_1, gen6_2, gen6_3};
+	funcs[7] = {gen7_1, gen7_2, gen7_3};
+	funcs[8] = {gen8_1, gen8_2};
+	size_t lastFuncMod = 9;
+
+	FuncOut out;
+	int nodes = 0;
+	while (nodes < size) {
+		size_t pos = std::random_devices::semirandom() % lastFuncMod;
+		if (funcs.find(pos) == funcs.end())
+			continue;
+		mergeGraphs(out, funcs.at(pos).at(std::random_devices::semirandom() % funcs.at(pos).size())());
+		nodes += pos;
+	}
+	return out.graph;
+}
+
+UndirectedGraph RandomGraphFactory::generateUndirectedIsomorphicGraph(const UndirectedGraph &graph)
+{
+	UndirectedGraph out;
+	std::unordered_map<Node, Node> mapping;
+
+	for (const Node &n : graph.getNodes()) {
+		Node node(std::string(n.getName()) + "i");
+		mapping.insert({n, node});
+		out.addNode(node);
+	}
+
+	for (const UndirectedEdge &e : graph.getEdges()) {
+		UndirectedEdge edge(mapping.at(e.getFirstNode()), mapping.at(e.getSecondNode()), e.getName());
+		out.addEdge(edge);
+	}
+
+	const std::set<UndirectedEdge> &edges = out.getEdges();
+
+	for (size_t i = 0; i < edges.size(); ++i) {
+		size_t pos = std::random_devices::semirandom() % edges.size();
+		UndirectedEdge e = *std::next(edges.begin(), pos);
+		out.removeEdge(e);
+		if (std::random_devices::semirandom() % 2)
+			e = UndirectedEdge(e.getSecondNode(), e.getFirstNode(), e.getName());
+		out.addEdge(e);
+	}
+
+	return out;
+}
+
+} // namespace generate
+
+} // namespace graph
diff --git a/alib2algo/src/graph/generate/RandomGraphFactory.h b/alib2algo/src/graph/generate/RandomGraphFactory.h
new file mode 100644
index 0000000000000000000000000000000000000000..7727f1fa9fe2d1a0a301c48a3bc5ae9cdded273a
--- /dev/null
+++ b/alib2algo/src/graph/generate/RandomGraphFactory.h
@@ -0,0 +1,25 @@
+#ifndef RANDOM_GRAPH_FACTORY_H_
+#define RANDOM_GRAPH_FACTORY_H_
+
+#include <graph/directed/DirectedGraph.h>
+#include <graph/undirected/UndirectedGraph.h>
+
+namespace graph
+{
+
+namespace generate
+{
+
+class RandomGraphFactory
+{
+public:
+	// size - rough number of nodes, may be more but not less
+	static UndirectedGraph generateUndirectedTriConnectedPlanarGraph(int size);
+	static UndirectedGraph generateUndirectedIsomorphicGraph(const UndirectedGraph &graph);
+};
+
+} // namespace generate
+
+} // namespace grammar
+
+#endif // RANDOM_GRAPH_FACTORY_H_
diff --git a/alib2algo/src/graph/isomorphism/HopcroftDebug.h b/alib2algo/src/graph/isomorphism/HopcroftDebug.h
new file mode 100644
index 0000000000000000000000000000000000000000..0863c9ace7af2e6d1ef6297f100e99683fbb31e5
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftDebug.h
@@ -0,0 +1,13 @@
+#ifndef GRAPH_HOPCROFT_DEBUG_H_
+#define GRAPH_HOPCROFT_DEBUG_H_
+
+#ifndef NDEBUG
+#include <cassert>
+#define H_ASSERT(x) assert(x)
+#define H_DEBUG(x) std::cout << std::boolalpha << x << std::endl
+#else
+#define H_ASSERT(x)
+#define HOPCROFT_DEBUG(str)
+#endif
+
+#endif // GRAPH_HOPCROFT_DEBUG_H_
diff --git a/alib2algo/src/graph/isomorphism/HopcroftGraph.cpp b/alib2algo/src/graph/isomorphism/HopcroftGraph.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..48792b990135cc5e9fac85cf4d59f2ab2ad11852
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftGraph.cpp
@@ -0,0 +1,806 @@
+#include "HopcroftGraph.h"
+#include "HopcroftDebug.h"
+
+#include <algorithm>
+
+#include <exception/AlibException.h>
+
+#include <iostream>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+namespace hopcroft
+{
+
+Graph::Graph(const UndirectedGraph &graph)
+{
+	init(graph);
+}
+
+Graph::~Graph()
+{
+	for (Face *f : faces)
+		delete f;
+
+	for (Vertex *v : vertices)
+		delete v;
+
+	for (EdgeEnd *e : edge_ends)
+		delete e;
+}
+
+// In: e1 <cw> [newedge] <cw> e2->otheredge <cw> e1
+// Out: newedge
+EdgeEnd *Graph::addInnerEdge(EdgeEnd *e1, EdgeEnd *e2)
+{
+	H_ASSERT(e1->vertex != e2->vertex);
+	H_ASSERT(e1->otheredge->ccwedge == e2);
+	H_ASSERT(std::find(edge_ends.begin(), edge_ends.end(), e1) != edge_ends.end());
+	H_ASSERT(std::find(edge_ends.begin(), edge_ends.end(), e2) != edge_ends.end());
+
+	//H_DEBUG("add inner edge" << e1->dump() << e2->dump());
+
+	EdgeEnd *enew = new EdgeEnd;
+	enew->vertex = e2->vertex;
+
+	EdgeEnd *enewother = new EdgeEnd;
+	enewother->vertex = e1->vertex;
+
+	enew->otheredge = enewother;
+	enewother->otheredge = enew;
+
+	enew->cwedge = e2->otheredge;
+	enew->ccwedge = e2->ccwedge;
+	e2->ccwedge->otheredge->cwedge = enewother;
+	e2->ccwedge = enewother;
+	if (e2->cwedge == e2->otheredge)
+		e2->cwedge = enewother;
+
+	enewother->cwedge = e1->cwedge;
+	enewother->ccwedge = e1->otheredge;
+	e1->cwedge->otheredge->ccwedge = enew;
+	e1->cwedge = enew;
+	if (e1->ccwedge == e1->otheredge)
+		e1->ccwedge = enew;
+
+	if (!enew->isMultiEdge()) {
+		enew->vertex->degree++;
+		enewother->vertex->degree++;
+	}
+
+	Face *facenew = new Face;
+	facenew->edges = 3;
+	facenew->fedge = e2;
+
+	e1->otheredge->ccwface = facenew;
+	enewother->ccwface = facenew;
+	e2->ccwface = facenew;
+
+	Face *faceold = enew->ccwedge->ccwface;
+	faceold->edges--;
+	faceold->fedge = enew;
+	enew->ccwface = faceold;
+
+	faces.push_back(facenew);
+	edge_ends.push_back(enew);
+	edge_ends.push_back(enewother);
+
+	return enew;
+}
+
+// In: e distinguished edge of 3-face
+// Out: face edges are reconnected to the new vertex
+void Graph::addAndReconnectVertexIntoTriangularFace(EdgeEnd *e)
+{
+	Vertex *v = new Vertex;
+	EdgeEnd *f = e->ccwedge;
+	EdgeEnd *g = f->ccwedge;
+	Face *oldface = e->ccwface;
+
+	H_ASSERT(e == g->ccwedge);
+
+	e->cwedge->otheredge->ccwedge = f;
+	f->cwedge->otheredge->ccwedge = g;
+	g->cwedge->otheredge->ccwedge = e;
+
+	e->cwedge = f->otheredge;
+	e->ccwedge = g->otheredge;
+	f->cwedge = g->otheredge;
+	f->ccwedge = e->otheredge;
+	g->cwedge = e->otheredge;
+	g->ccwedge = f->otheredge;
+
+	v->degree = 3;
+	v->vedge = e->otheredge;
+
+	e->otheredge->vertex->updateDegree();
+	f->otheredge->vertex->updateDegree();
+	g->otheredge->vertex->updateDegree();
+
+	e->ccwface = g->otheredge->ccwface;
+	f->ccwface = e->otheredge->ccwface;
+	g->ccwface = f->otheredge->ccwface;
+
+	e->ccwface->edges++;
+	f->ccwface->edges++;
+	g->ccwface->edges++;
+
+	e->ccwface->fedge = e;
+	f->ccwface->fedge = f;
+	g->ccwface->fedge = g;
+
+	removeFace(oldface);
+	vertices.push_back(v);
+}
+
+// In: v -> w edge-end (v = the vertex)
+// Out: w remains
+void Graph::removeDistinguishedEdgeAndVertex(EdgeEnd *e)
+{
+	H_ASSERT(std::find(edge_ends.begin(), edge_ends.end(), e) != edge_ends.end());
+
+	//H_DEBUG("remove distinguished edge and vertex" << e->dump());
+
+	EdgeEnd *o = e->otheredge;
+
+	Vertex *w = e->vertex;
+	Vertex *v = o->vertex;
+
+	EdgeEnd *wl = e->cwedge;
+	EdgeEnd *wk = e->ccwedge;
+	EdgeEnd *vl = o->cwedge;
+	EdgeEnd *vk = o->ccwedge;
+
+	if (wk == o)
+		wk = vk;
+
+	if (wl == o)
+		wl = vl;
+
+	v->forEachEdge([=](EdgeEnd *e) {
+		e->otheredge->vertex = w;
+		return true;
+	});
+
+	wl->otheredge->ccwedge = vk;
+	wk->otheredge->cwedge = vl;
+	vl->otheredge->ccwedge = wk;
+	vk->otheredge->cwedge = wl;
+
+	w->vedge = wk;
+	w->updateDegree();
+
+	w->forEachEdge([=](EdgeEnd *e) {
+		e->vertex->updateDegree();
+		return true;
+	});
+
+	wk->ccwface->edges--;
+	wk->ccwface->fedge = wk;
+	vk->ccwface->edges--;
+	vk->ccwface->fedge = vk;
+
+	removeVertex(v);
+	removeEdgeEnd(e);
+	removeEdgeEnd(o);
+}
+
+// In: v - exception vertex, edge e points to 4 degree vertex
+void Graph::removeFourDegreeExceptionVertex(Vertex *v, EdgeEnd *e)
+{
+	H_ASSERT(v->degree == 4);
+	H_ASSERT(e->vertex->degree == 4);
+	H_ASSERT(!e->isMultiEdge());
+
+	//H_DEBUG("remove four degree exception vertex" << v->dump() << e->dump());
+
+	// Out: e <ccw> dest
+	auto reconnectEdge = [=](EdgeEnd *e, EdgeEnd *dest) {
+		e->vertex = dest->otheredge->vertex;
+		e->ccwedge = dest;
+		e->cwedge = dest->otheredge->cwedge;
+		dest->otheredge->cwedge = e->otheredge;
+		if (dest->otheredge->ccwedge == dest->otheredge)
+			dest->otheredge->ccwedge = e->otheredge;
+		e->cwedge->otheredge->ccwedge = e->otheredge;
+	};
+
+	EdgeEnd *er = e;
+	EdgeEnd *es = er->otheredge->cwedge;
+	EdgeEnd *et = es->otheredge->cwedge;
+	EdgeEnd *eu = et->otheredge->cwedge;
+
+	EdgeEnd *destroyed1 = es->otheredge;
+	EdgeEnd *destroyed2 = eu->otheredge;
+
+	reconnectEdge(et->otheredge, eu->otheredge);
+	reconnectEdge(er->otheredge, es->otheredge);
+
+	// Merge es + eu
+	es->otheredge = eu;
+	eu->otheredge = es;
+	es->vertex->vedge = eu;
+	eu->vertex->vedge = es;
+	es->cwedge->otheredge->ccwedge = eu;
+	es->ccwedge->otheredge->cwedge = eu;
+	eu->cwedge->otheredge->ccwedge = es;
+	eu->ccwedge->otheredge->cwedge = es;
+
+	// Update faces
+	Face *face1 = er->otheredge->ccwface;
+	Face *face2 = et->otheredge->ccwface;
+	Face *face3 = er->ccwface;
+	Face *face4 = et->ccwface;
+
+	eu->ccwface = face1;
+	es->ccwface = face2;
+	face1->fedge = eu;
+	face2->fedge = es;
+	face3->fedge = er;
+	face4->fedge = et;
+	face3->edges--;
+	face4->edges--;
+
+	// Update degrees
+	er->vertex->updateDegree();
+	es->vertex->updateDegree();
+	et->vertex->updateDegree();
+	eu->vertex->updateDegree();
+
+	removeVertex(v);
+	removeEdgeEnd(destroyed1);
+	removeEdgeEnd(destroyed2);
+}
+
+void Graph::removeEdge(EdgeEnd *e, bool removeVertices)
+{
+	H_ASSERT(std::find(edge_ends.begin(), edge_ends.end(), e) != edge_ends.end());
+
+	//H_DEBUG("remove edge" << e->dump() << e->otheredge->dump());
+
+	EdgeEnd *o = e->otheredge;
+
+	if (!e->isMultiEdge()) {
+		e->vertex->degree--;
+		// In case of loops
+		if (e->vertex != o->vertex)
+			o->vertex->degree--;
+	}
+
+	// Update edges
+	e->cwedge->otheredge->ccwedge = e->ccwedge;
+	if (e->cwedge->otheredge->cwedge == o)
+		e->cwedge->otheredge->cwedge = e->cwedge->otheredge;
+
+	e->ccwedge->otheredge->cwedge = e->cwedge;
+	if (e->ccwedge->otheredge->ccwedge == o)
+		e->ccwedge->otheredge->ccwedge = e->ccwedge->otheredge;
+
+	o->cwedge->otheredge->ccwedge = o->ccwedge;
+	if (o->cwedge->otheredge->cwedge == e)
+		o->cwedge->otheredge->cwedge = o->cwedge->otheredge;
+
+	o->ccwedge->otheredge->cwedge = o->cwedge;
+	if (o->ccwedge->otheredge->ccwedge == e)
+		o->ccwedge->otheredge->ccwedge = o->ccwedge->otheredge;
+
+	// Update vedges
+	if (e->vertex->vedge == o)
+		e->vertex->vedge = e->cwedge;
+
+	if (o->vertex->vedge == e)
+		o->vertex->vedge = o->cwedge;
+
+	// Handle zero degree vertices
+	if (removeVertices) {
+		if (e->vertex->degree == 0)
+			removeVertex(e->vertex);
+
+		if (o->vertex->degree == 0)
+			removeVertex(o->vertex);
+	} else {
+		if (e->vertex->degree == 0)
+			e->vertex->vedge = nullptr;
+
+		if (o->vertex->degree == 0)
+			o->vertex->vedge = nullptr;
+	}
+
+	// Update faces
+	Face *face1 = e->ccwface;
+	Face *face2 = o->ccwface;
+
+	if (face1 != face2) {
+		face1->edges = 0;
+		face1->fedge = (e->ccwedge == e || e->ccwedge == o) ? o->ccwedge : e->ccwedge;
+		EdgeEnd *edge = face1->fedge;
+		do {
+			face1->edges++;
+			edge->ccwface = face1;
+			edge = edge->ccwedge;
+		} while (edge != face1->fedge);
+		removeFace(face2);
+	} else {
+		face1->edges -= 2;
+		face1->fedge = o->ccwedge == o ? e->ccwedge : o->ccwedge;
+		if (face1->edges == 0)
+			removeFace(face1);
+	}
+
+	// Finally remove edge
+	removeEdgeEnd(e);
+	removeEdgeEnd(o);
+}
+
+void Graph::removeFace(Face *f)
+{
+	H_ASSERT(std::find(faces.begin(), faces.end(), f) != faces.end());
+
+	faces.erase(std::find(faces.begin(), faces.end(), f));
+	delete f;
+}
+
+void Graph::removeVertex(Vertex *v)
+{
+	H_ASSERT(std::find(vertices.begin(), vertices.end(), v) != vertices.end());
+
+	vertices.erase(std::find(vertices.begin(), vertices.end(), v));
+	delete v;
+}
+
+void Graph::removeEdgeEnd(EdgeEnd *e)
+{
+	H_ASSERT(std::find(edge_ends.begin(), edge_ends.end(), e) != edge_ends.end());
+
+	edge_ends.erase(std::find(edge_ends.begin(), edge_ends.end(), e));
+	delete e;
+}
+
+std::vector<std::vector<EdgeEnd*>> Graph::findClumps(Vertex *v, Vertex *w) const
+{
+	std::vector<EdgeEnd*> clump;
+	std::vector<std::vector<EdgeEnd*>> clumps;
+
+	v->forEachEdge([&](EdgeEnd *e) {
+		if (e->vertex == w && e->otheredge->vertex == v && (clump.empty() || clump.back()->cwedge == e->otheredge)) {
+			clump.push_back(e);
+		} else if (!clump.empty()) {
+			clumps.push_back(clump);
+			clump.clear();
+			clump.push_back(e);
+		}
+		return true;
+	});
+
+	if (!clump.empty())
+		clumps.push_back(clump);
+
+	if (clumps.empty())
+		return clumps;
+
+	// Merge potential <-clump1) .. (clump2) .. (clump1->
+	if (clumps.size() > 1 && clumps.back().back()->cwedge == clumps.front().front()->otheredge) {
+		clumps[0].insert(clumps[0].end(), clumps.back().begin(), clumps.back().end());
+		clumps.pop_back();
+	}
+
+	// Remove size-one clumps
+	for (size_t i = 0; i < clumps.size(); ++i)
+		if (clumps[i].size() < 2)
+			clumps.erase(clumps.begin() + i--);
+
+	return clumps;
+}
+
+void Graph::init(const UndirectedGraph &graph)
+{
+	std::unordered_map<Node, Vertex*> v_map;
+	std::unordered_map<UndirectedEdge, EdgeEnd*> e_map;
+
+	for (const Node &node : graph.getNodes()) {
+		Vertex *v = new Vertex;
+		v->out = std::string(node.getName());
+		vertices.push_back(v);
+		v_map[node] = v;
+	}
+
+	for (const UndirectedEdge &edge : graph.getEdges()) {
+		EdgeEnd *e1 = new EdgeEnd;
+		e1->vertex = v_map[edge.getSecondNode()];
+		edge_ends.push_back(e1);
+
+		EdgeEnd *e2 = new EdgeEnd;
+		e2->vertex = v_map[edge.getFirstNode()];
+		edge_ends.push_back(e2);
+
+		e1->otheredge = e2;
+		e2->otheredge = e1;
+		e1->vertex->vedge = e2;
+		e2->vertex->vedge = e1;
+
+		e_map[edge] = e1;
+	}
+
+	// Embedding
+	std::unordered_map<Node, std::vector<EdgeEnd*>> sortedNodeEdges;
+
+	for (const auto &p : graph.getEmbedding()) {
+		Vertex *vertex = v_map[p.first];
+		std::vector<EdgeEnd*> sortedEdges;
+		std::set<UndirectedEdge> neighborEdges = graph.neighborEdges(p.first);
+
+		// Multi-edges must not cross the vector boundary in embedding
+		// eg. (n1, n2, n3, n4, n1) must be changed to (n1, n1, n2, n3, n4)
+		auto reorder = [&](const std::vector<Node> &v) {
+			if (v_map.size() == 1 || v.front() != v.back())
+				return v;
+
+			std::vector<Node> out = v;
+			while (out.front() == out.back()) {
+				out.insert(out.begin(), out.back());
+				out.pop_back();
+			}
+			return out;
+		};
+
+		for (const Node &node : reorder(p.second)) {
+			Vertex *v = v_map[node];
+			EdgeEnd *edge = nullptr;
+
+			for (const UndirectedEdge &edg : neighborEdges) {
+				EdgeEnd *e = e_map[edg];
+				if (e->vertex == vertex)
+					e = e->otheredge;
+				if (e->vertex == v) {
+					edge = e;
+					const Node &othernode = edg.getFirstNode() == p.first ? edg.getSecondNode() : edg.getFirstNode();
+					auto &vec = sortedNodeEdges[othernode];
+					for (size_t i = vec.size(); i-- > 0; ) {
+						bool isMultiEdge = false;
+						EdgeEnd *e = vec[i]->otheredge;
+						if (e->ccwedge && e->ccwedge->vertex == e->otheredge->vertex)
+							isMultiEdge = true;
+						if (e->cwedge && e->cwedge->vertex == e->otheredge->vertex)
+							isMultiEdge = true;
+						if (isMultiEdge && e->vertex == v && e->otheredge->vertex == vertex) {
+							vec.erase(vec.begin() + i);
+							edge = e;
+							break;
+						}
+					}
+
+					neighborEdges.erase(edg);
+					break;
+				}
+			}
+
+			H_ASSERT(edge);
+			sortedEdges.push_back(edge);
+		}
+
+		sortedNodeEdges[p.first] = sortedEdges;
+
+		for (size_t i = 0; i < sortedEdges.size(); ++i) {
+			EdgeEnd *e1 = sortedEdges[i];
+			EdgeEnd *e2 = i == sortedEdges.size() - 1 ? sortedEdges[0] : sortedEdges[i + 1];
+
+			H_ASSERT(e1 != e2);
+			H_ASSERT(!e1->isLoopEdge());
+
+			// CW                            // CCW
+			e1->otheredge->ccwedge = e2;     // e1->otheredge->cwedge = e2;
+			e2->otheredge->cwedge = e1;      // e2->otheredge->ccwedge = e1;
+		}
+	}
+
+	if (sortedNodeEdges.empty() && !e_map.empty())
+		throw exception::AlibException("HopcroftWong: Graph is not embedded");
+
+	for (Vertex *vertex : vertices) {
+		vertex->degree = 1;
+		vertex->updateDegree();
+	}
+
+	// Init faces
+	for (EdgeEnd *e : edge_ends) {
+		if (e->ccwface)
+			continue;
+		Face *f = new Face;
+		faces.push_back(f);
+		f->fedge = e;
+		EdgeEnd *next = e;
+		do {
+			f->edges++;
+			next->ccwface = f;
+			next = next->ccwedge;
+		} while (next != e);
+	}
+}
+
+void Graph::dump()
+{
+	std::cout << std::endl;
+	std::cout << "Vertices: " << vertices.size();
+	std::cout << " Edges: " << edge_ends.size() / 2;
+	std::cout << std::endl;
+	for (Vertex *v : vertices) {
+		std::cout << v->dump() << "(" << v->degree << "):";
+		v->forEachEdge([this](EdgeEnd *e) {
+			std::cout << e->dump() << " (" << std::distance(edge_ends.begin(), std::find(edge_ends.begin(), edge_ends.end(), e)) << ")";
+			return true;
+		});
+		std::cout << std::endl;
+	}
+	std::cout << std::endl;
+}
+
+void Graph::dumpLabels()
+{
+	std::cout << std::endl;
+	std::cout << "Vertices:" << std::endl;
+	for (Vertex *v : vertices)
+		std::cout << " " << v->vlabel;
+	std::cout << std::endl;
+
+	std::cout << "Edges:" << std::endl;
+	for (Vertex *v : vertices) {
+		std::cout << v->dump() << ":";
+		v->forEachEdge([this](EdgeEnd *e) {
+			std::cout << " " << e->elabel << "-" << e->otheredge->elabel;
+			return true;
+		});
+		std::cout << std::endl;
+	}
+	std::cout << std::endl;
+}
+
+Vertex::Vertex()
+	: vlabel(1)
+	, degree(0)
+	, vedge(nullptr)
+{
+}
+
+bool Vertex::isKnot() const
+{
+	if (degree != 1)
+		return false;
+
+	int loops = 0;
+	bool res = true;
+
+	forEachEdge([&](EdgeEnd *e) {
+		if (!e->isLoopEdge()) {
+			res = false;
+			return false;
+		}
+		loops++;
+		return true;
+	});
+
+	return res && loops == 1;
+}
+
+bool Vertex::isCorolla() const
+{
+	if (degree != 1)
+		return false;
+
+	int loops = 0;
+	bool res = true;
+
+	forEachEdge([&](EdgeEnd *e) {
+		if (!e->isLoopEdge()) {
+			res = false;
+			return false;
+		}
+		loops++;
+		return true;
+	});
+
+	return res && loops > 1;
+}
+
+bool Vertex::isStar() const
+{
+	int edges = 0;
+	bool res = true;
+
+	forEachEdge([&](EdgeEnd *e) {
+		if (!e->isSpoke())
+			res = false;
+		edges++;
+		return res;
+	});
+
+	return res && edges > 1;
+}
+
+bool Vertex::isDumbell() const
+{
+	int edges = 0;
+	bool res = true;
+
+	forEachEdge([&](EdgeEnd *e) {
+		if (!e->isSpoke())
+			res = false;
+		edges++;
+		return res;
+	});
+
+	return res && edges == 1;
+}
+
+bool Vertex::isLoopVertex() const
+{
+	bool res = false;
+
+	forEachEdge([&](EdgeEnd *e) {
+		res = e->isLoopEdge();
+		return !res;
+	});
+
+	return res;
+}
+
+bool Vertex::isSpokeCenter() const
+{
+	bool res = false;
+
+	forEachEdge([&](EdgeEnd *e) {
+		res = e->isSpoke();
+		return !res;
+	});
+
+	return res;
+}
+
+bool Vertex::hasNonLoopEdge() const
+{
+	bool res = false;
+
+	forEachEdge([&](EdgeEnd *e) {
+		res = !e->isLoopEdge();
+		return !res;
+	});
+
+	return res;
+}
+
+bool Vertex::hasNonSpokeEdge() const
+{
+	bool res = false;
+
+	forEachEdge([&](EdgeEnd *e) {
+		res = !e->isSpoke();
+		return !res;
+	});
+
+	return res;
+}
+
+bool Vertex::isDegreeIsolated() const
+{
+	bool res = true;
+
+	forEachEdge([&](EdgeEnd *e) {
+		if (e->vertex->degree == degree)
+			res = false;
+		return res;
+	});
+
+	return res;
+}
+
+std::vector<Face*> Vertex::neighborFaces(Direction direction) const
+{
+	std::vector<Face*> out;
+
+	forEachEdge([&](EdgeEnd *e) {
+		out.push_back(e->ccwface);
+		return true;
+	}, direction);
+
+	return out;
+}
+
+std::vector<EdgeEnd*> Vertex::neighborEdges(Direction direction) const
+{
+	std::vector<EdgeEnd*> out;
+
+	forEachEdge([&](EdgeEnd *e) {
+		out.push_back(e);
+		return true;
+	}, direction);
+
+	return out;
+}
+
+void Vertex::forEachEdge(std::function<bool(EdgeEnd*)> callback, Direction direction) const
+{
+	// Current edge can be removed from within the callback (without removing vertices)
+
+	if (!vedge)
+		return;
+
+	EdgeEnd *s = vedge;
+	EdgeEnd *e = vedge;
+	do {
+		EdgeEnd *n = direction == ClockWise ? e->otheredge->ccwedge : e->otheredge->cwedge;
+		if (!callback(e))
+			return;
+		e = n;
+	} while (e != s && degree > 0);
+}
+
+void Vertex::updateDegree()
+{
+	std::set<Vertex*> nghbrs;
+
+	forEachEdge([&](EdgeEnd *e) {
+		nghbrs.insert(e->vertex);
+		return true;
+	});
+
+	degree = nghbrs.size();
+}
+
+std::string Vertex::dump() const
+{
+	return " " + out;
+}
+
+EdgeEnd::EdgeEnd()
+	: elabel(2)
+	, vertex(nullptr)
+	, ccwface(nullptr)
+	, ccwedge(nullptr)
+	, cwedge(nullptr)
+	, otheredge(nullptr)
+{
+}
+
+bool EdgeEnd::isSpoke() const
+{
+	return !isMultiEdge() && (vertex->degree == 1 || otheredge->vertex->degree == 1);
+}
+
+bool EdgeEnd::isLoopEdge() const
+{
+	return vertex == otheredge->vertex;
+}
+
+bool EdgeEnd::isMultiEdge() const
+{
+	bool res = false;
+
+	otheredge->vertex->forEachEdge([&](EdgeEnd *e) {
+		if (e == this || e == otheredge)
+			return true;
+		res = e->vertex == vertex && e->otheredge->vertex == otheredge->vertex;
+		return !res;
+	});
+
+	return res;
+}
+
+std::string EdgeEnd::dump() const
+{
+	return " " + otheredge->vertex->out + "->" + vertex->out;
+}
+
+Face::Face()
+	: edges(0)
+	, fedge(nullptr)
+{
+}
+
+} // namespace hopcroft
+
+} // namespace isomorphism
+
+} // namespace graph
diff --git a/alib2algo/src/graph/isomorphism/HopcroftGraph.h b/alib2algo/src/graph/isomorphism/HopcroftGraph.h
new file mode 100644
index 0000000000000000000000000000000000000000..5242c9d6069dc8c44080b6acec35aa44eca214ac
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftGraph.h
@@ -0,0 +1,129 @@
+#ifndef GRAPH_HOPCROFT_GRAPH_H_
+#define GRAPH_HOPCROFT_GRAPH_H_
+
+#include <vector>
+
+#include <graph/Graph.h>
+#include <graph/directed/DirectedGraph.h>
+#include <graph/undirected/UndirectedGraph.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+namespace hopcroft
+{
+
+class Vertex;
+class EdgeEnd;
+class Face;
+
+enum Direction {
+	ClockWise,
+	CounterClockWise
+};
+
+class Graph
+{
+public:
+	explicit Graph(const UndirectedGraph &graph);
+	~Graph();
+
+	// Reduction-specific modifications
+	EdgeEnd *addInnerEdge(EdgeEnd *e1, EdgeEnd *e2);
+	void removeDistinguishedEdgeAndVertex(EdgeEnd *e);
+	void addAndReconnectVertexIntoTriangularFace(EdgeEnd *e);
+	void removeFourDegreeExceptionVertex(Vertex *v, EdgeEnd *e);
+
+	void removeEdge(EdgeEnd *e, bool removeVertices = true);
+
+	void removeFace(Face *f);
+	void removeVertex(Vertex *v);
+	void removeEdgeEnd(EdgeEnd *e);
+
+	std::vector<std::vector<EdgeEnd*>> findClumps(Vertex *v, Vertex *w) const;
+
+	void dump();
+	void dumpLabels();
+
+	std::vector<Face*> faces;
+	std::vector<Vertex*> vertices;
+	std::vector<EdgeEnd*> edge_ends;
+
+private:
+	Graph(const Graph &other) = delete;
+
+	void init(const UndirectedGraph &graph);
+};
+
+class Vertex
+{
+public:
+	explicit Vertex();
+
+	bool isKnot() const;
+	bool isCorolla() const;
+
+	bool isStar() const;
+	bool isDumbell() const;
+	bool isLoopVertex() const;
+	bool isSpokeCenter() const;
+
+	bool hasNonLoopEdge() const;
+	bool hasNonSpokeEdge() const;
+	bool isDegreeIsolated() const;
+
+	std::vector<Face*> neighborFaces(Direction d = ClockWise) const;
+	std::vector<EdgeEnd*> neighborEdges(Direction d = ClockWise) const;
+	void forEachEdge(std::function<bool(EdgeEnd*)> cb, Direction d = ClockWise) const;
+
+	void updateDegree();
+
+	int vlabel;
+	int degree;
+	EdgeEnd *vedge;
+
+	// Debug
+	std::string dump() const;
+	std::string out;
+};
+
+class EdgeEnd
+{
+public:
+	explicit EdgeEnd();
+
+	bool isSpoke() const;
+	bool isLoopEdge() const;
+
+	bool isMultiEdge() const;
+
+	int elabel;
+	Vertex *vertex;
+	Face *ccwface;
+	EdgeEnd *ccwedge;
+	EdgeEnd *cwedge;
+	EdgeEnd *otheredge; // other edge-end
+
+	// Debug
+	std::string dump() const;
+};
+
+class Face
+{
+public:
+	explicit Face();
+
+	int edges;
+	EdgeEnd *fedge;
+};
+
+} // namespace hopcroft
+
+} // namespace isomorphism
+
+} // namespace graph
+
+#endif // GRAPH_HOPCROFT_GRAPH_H_
diff --git a/alib2algo/src/graph/isomorphism/HopcroftImpl.cpp b/alib2algo/src/graph/isomorphism/HopcroftImpl.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bfb7ddfc226d09d0c78cc9386d2a68fc18b2ecbe
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftImpl.cpp
@@ -0,0 +1,1500 @@
+#include "HopcroftImpl.h"
+#include "HopcroftDebug.h"
+
+#include <tuple>
+#include <iostream>
+#include <algorithm>
+
+#include <exception/AlibException.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+namespace hopcroft
+{
+
+// Helpers
+
+#define ARG_REF(x) decltype(x##1) &x
+#define ARG_CREF(x) const decltype(x##1) &x
+
+#define REDUCTION_INIT \
+	int __flab = 0; \
+	(void)__flab;
+
+#define PREPARE_FUNCTION REDUCTION_INIT auto findReductions = [&]
+#define APPLY_FUNCTION auto applyReductions = [&]
+
+#define GET_LABELS(labels) \
+	( __flab++ == 0 \
+		? std::vector<int>(labels.begin(), labels.begin() + labels.size() / 2) \
+		: std::vector<int>(labels.begin() + labels.size() / 2 , labels.end()) \
+	)
+
+#define VECTOR_PARTITIONING(numbers, labels) \
+	labels = vectorPartitioning(numbers);
+
+#define CIRCLE_PARTITIONING(numbers, labels) \
+	labels = circlePartitioning(numbers);
+
+#define CHECK_REDUCTIONS(x) \
+	if (x##1 .size() != x##2 .size() || __ig1 != __ig2) { \
+		m_fail = true; \
+		return true; \
+	} \
+	if (x##1 .empty()) { \
+		return false; \
+	}
+
+#define CHECK_REDUCTIONS2(x1, x2) \
+	if (x1##1 .size() != x1##2 .size() || x2##1 .size() != x2##2 .size() || __ig1 != __ig2) { \
+		m_fail = true; \
+		return true; \
+	} \
+	if (x1##1 .empty() && x2##1 .empty()) { \
+		return false; \
+	}
+
+#define PREPARE_REDUCTIONS_X(x1, x2) \
+	saveInteger(); \
+	x1; \
+	int __ig1 = m_integer; \
+	restoreInteger(); \
+	x2; \
+	int __ig2 = m_integer; \
+
+#define PREPARE_REDUCTIONS2(x1, x2) \
+	PREPARE_REDUCTIONS_X( \
+		findReductions(x1##1, x2##1), \
+		findReductions(x1##2, x2##2))
+
+#define PREPARE_REDUCTIONS3(x1, x2, x3) \
+	PREPARE_REDUCTIONS_X( \
+		findReductions(x1##1, x2##1, x3##1), \
+		findReductions(x1##2, x2##2, x3##2))
+
+#define APPLY_REDUCTIONS_X(x1, x2) \
+	saveInteger(); \
+	x1; \
+	__ig1 = m_integer; \
+	restoreInteger(); \
+	x2; \
+	__ig2 = m_integer; \
+	m_fail = __ig1 != __ig2;
+
+#define APPLY_REDUCTIONS2(x1, x2) \
+	APPLY_REDUCTIONS_X( \
+		applyReductions(x1##1, x2##1), \
+		applyReductions(x1##2, x2##2))
+
+#define APPLY_REDUCTIONS3(x1, x2, x3) \
+	APPLY_REDUCTIONS_X( \
+		applyReductions(x1##1, x2##1, x3##1), \
+		applyReductions(x1##2, x2##2, x3##2))
+
+template<typename T, typename S>
+inline bool contains(const T &container, const S &value)
+{
+	return std::find(container.begin(), container.end(), value) != container.end();
+}
+
+template<typename T, typename S>
+inline int indexOf(const T &container, const S &value)
+{
+	auto it = std::find(container.begin(), container.end(), value);
+	if (it == container.end())
+		return -1;
+	return std::distance(container.begin(), it);
+}
+
+template<typename T>
+inline void deleteAll(const T &container)
+{
+	for (auto v : container)
+		delete v;
+}
+
+template<typename T>
+inline void deleteAllValues(const T &container)
+{
+	for (const auto &v : container)
+		delete v.second;
+}
+
+static const int SPECIAL_X = 12345;
+
+// HopcroftImpl
+
+HopcroftImpl::HopcroftImpl(Graph *g1, Graph *g2)
+	: m_g1(g1)
+	, m_g2(g2)
+	, m_fail(false)
+	, m_integer(3) // 1 - new vertices, 2 - new edges
+	, m_savedInteger(0)
+{
+}
+
+bool HopcroftImpl::testIsomorphism()
+{
+	while (doReductions()) {
+		if (m_g1->vertices.size() != m_g2->vertices.size())
+			m_fail = true;
+		if (m_g1->edge_ends.size() != m_g2->edge_ends.size())
+			m_fail = true;
+		if (m_g1->faces.size() != m_g2->faces.size())
+			m_fail = true;
+		if (m_fail)
+			break;
+	}
+
+	return !m_fail && testByDefinition();
+}
+
+int HopcroftImpl::newInteger()
+{
+	return m_integer++;
+}
+
+void HopcroftImpl::saveInteger()
+{
+	m_savedInteger = m_integer;
+}
+
+void HopcroftImpl::restoreInteger()
+{
+	m_integer = m_savedInteger;
+}
+
+// 2.A Vector partitioning
+std::vector<int> HopcroftImpl::vectorPartitioning(const std::vector<std::vector<int>> &tuples)
+{
+	std::vector<int> out(tuples.size(), -1);
+
+	for (size_t i = 0; i < out.size(); ++i) {
+		if (out[i] != -1)
+			continue;
+		int c = newInteger();
+		out[i] = c;
+		for (size_t j = 0; j < out.size(); ++j)
+			if (tuples[i] == tuples[j])
+				out[j] = c;
+	}
+
+	return out;
+}
+
+// 2.B Circle partitioning
+bool HopcroftImpl::circleIsomorphism(std::vector<int> a, std::vector<int> b)
+{
+	auto distinctIntegers = [](const std::vector<int> &circle) {
+		if (circle.empty())
+			return false;
+
+		int first = circle.at(0);
+		for (const int &val : circle)
+			if (val != first)
+				return true;
+
+		return false;
+	};
+
+	int savedInteger = m_integer;
+
+	while (distinctIntegers(a) || distinctIntegers(b)) {
+		std::vector<std::vector<int>> tuples;
+
+		auto constructTuples = [&tuples](const std::vector<int> &circle) {
+			for (size_t i = 0; i < circle.size(); ++i) {
+				size_t first = i;
+				size_t second = (i == circle.size() - 1) ? 0 : i + 1;
+				if (circle.at(first) != circle.at(second))
+					tuples.push_back({circle.at(first), circle.at(second)});
+			}
+		};
+
+		auto replaceComponents = [&tuples](std::vector<int> &circle, int idx, int cls) {
+			const std::vector<int> &comp = tuples.at(idx);
+			H_ASSERT(comp.size() == 2);
+
+			while (true) {
+				size_t first;
+				size_t second;
+				bool found = false;
+				for (size_t i = 0; i < circle.size(); ++i) {
+					first = i;
+					second = (i == circle.size() - 1) ? 0 : i + 1;
+					if (circle.at(first) == comp.at(0) && circle.at(second) == comp.at(1)) {
+						found = true;
+						break;
+					}
+				}
+				if (!found)
+					break;
+
+				circle[first] = cls;
+				circle.erase(circle.begin() + second);
+			}
+		};
+
+		constructTuples(a);
+		constructTuples(b);
+
+		const std::vector<int> &part = vectorPartitioning(tuples);
+		for (size_t i = 0; i < part.size(); ++i) {
+			int cls = part.at(i);
+			for (size_t j = i; j < part.size(); ++j) {
+				if (part.at(j) != cls)
+					continue;
+
+				replaceComponents(a, j, cls);
+				replaceComponents(b, j, cls);
+			}
+		}
+	}
+
+	m_integer = savedInteger;
+
+	return a == b;
+}
+
+std::vector<int> HopcroftImpl::circlePartitioning(const std::vector<std::vector<int>> &circles)
+{
+	std::vector<int> out(circles.size(), -1);
+
+	for (size_t i = 0; i < out.size(); ++i) {
+		if (out[i] != -1)
+			continue;
+		int c = newInteger();
+		out[i] = c;
+		for (size_t j = 0; j < out.size(); ++j)
+			if (circleIsomorphism(circles[i], circles[j]))
+				out[j] = c;
+	}
+
+	return out;
+}
+
+bool HopcroftImpl::doReductions()
+{
+	// Reductions returns true when a reduction is applied or set flag when fail
+	// -> loop until there is no available reduction
+	if (reductionA() || m_fail)
+		return true;
+	if (reductionB() || m_fail)
+		return true;
+	if (reductionC() || m_fail)
+		return true;
+	return false;
+}
+
+bool HopcroftImpl::reductionA()
+{
+	// Removal of loops and degree one vertices
+	if (removeLoops())
+		return true;
+	if (removeCorollas())
+		return true;
+	if (removeKnots())
+		return true;
+	if (removeDegreeOneVertices())
+		return true;
+	return false;
+}
+
+// A(i) Removing loops from loop vertices with non-loop edges
+bool HopcroftImpl::removeLoops()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 3-tuple
+	std::vector<EdgeEnd*> loops_g1;
+	std::vector<EdgeEnd*> loops_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(loops_g)) {
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isLoopVertex() && vertex->hasNonLoopEdge()) {
+				vertex->forEachEdge([&](EdgeEnd *edge) {
+					if (edge->isLoopEdge() && (edge->cwedge == edge)) {
+						EdgeEnd *f = edge->otheredge;
+						EdgeEnd *e = f->cwedge;
+						loops_g.push_back(f);
+						numbers.push_back({e->otheredge->elabel, f->elabel, f->otheredge->elabel});
+					}
+					return true;
+				});
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(loops_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			loops_g[i]->cwedge->otheredge->elabel = labels_g[i];
+
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			graph->removeEdge(loops_g[i]);
+	};
+
+	PREPARE_REDUCTIONS2(m_g, loops_g);
+	CHECK_REDUCTIONS(loops_g);
+	H_DEBUG("Apply removeLoops");
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, loops_g);
+	return true;
+}
+
+// A(ii) Removing corollas
+bool HopcroftImpl::removeCorollas()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<Vertex*> corollas_g1;
+	std::vector<Vertex*> corollas_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(corollas_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isCorolla()) {
+				numbers.resize(numbers.size() + 1);
+				vertex->forEachEdge([&](EdgeEnd *e) {
+					EdgeEnd *f = e->cwedge;
+					if (f != e->otheredge) {
+						e = e->otheredge;
+						f = e->cwedge;
+					}
+					numbers[ni].push_back(vertex->vlabel);
+					numbers[ni].push_back(e->elabel);
+					numbers[ni].push_back(f->elabel);
+					return true;
+				});
+				ni++;
+				corollas_g.push_back(vertex);
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(corollas_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			Vertex *corolla = corollas_g[i];
+			while (corolla->vedge)
+				graph->removeEdge(corolla->vedge, /* removeVertices */ false);
+			corolla->vlabel = labels_g[i];
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, corollas_g);
+	CHECK_REDUCTIONS(corollas_g);
+	H_DEBUG("Apply removeCorollas");
+	CIRCLE_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, corollas_g);
+	return true;
+}
+
+// A(iii) Removing knots
+bool HopcroftImpl::removeKnots()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<Vertex*> knots_g1;
+	std::vector<Vertex*> knots_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(knots_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isKnot()) {
+				numbers.resize(numbers.size() + 1);
+
+				numbers[ni].push_back(vertex->vlabel);
+				numbers[ni].push_back(vertex->vedge->elabel);
+				numbers[ni].push_back(vertex->vedge->otheredge->elabel);
+				numbers[ni].push_back(vertex->vlabel);
+				numbers[ni].push_back(vertex->vedge->otheredge->elabel);
+				numbers[ni].push_back(vertex->vedge->elabel);
+
+				ni++;
+				knots_g.push_back(vertex);
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(knots_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			Vertex *knot = knots_g[i];
+			graph->removeEdge(knot->vedge, /* removeVertices */ false);
+			knot->vlabel = labels_g[i];
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, knots_g);
+	CHECK_REDUCTIONS(knots_g);
+	H_DEBUG("Apply removeKnots");
+	CIRCLE_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, knots_g);
+	return true;
+}
+
+// A(iv) Removing degree one vertices
+bool HopcroftImpl::removeDegreeOneVertices()
+{
+	if (removeSpokeCenters() || m_fail)
+		return true;
+	if (removeStars() || m_fail)
+		return true;
+	if (removeDumbells() || m_fail)
+		return true;
+	return false;
+}
+
+bool HopcroftImpl::removeSpokeCenters()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 3-tuple
+	std::vector<EdgeEnd*> spokes_g1;
+	std::vector<EdgeEnd*> spokes_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(spokes_g)) {
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isSpokeCenter() && vertex->hasNonSpokeEdge()) {
+				vertex->forEachEdge([&](EdgeEnd *edge) {
+					if (edge->isSpoke()) {
+						EdgeEnd *f = edge->otheredge;
+						EdgeEnd *e = f->cwedge;
+						spokes_g.push_back(f);
+						numbers.push_back({e->otheredge->elabel, f->elabel, f->otheredge->elabel});
+						H_ASSERT(!contains(spokes_g, f->otheredge));
+					}
+					return true;
+				});
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(spokes_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			spokes_g[i]->cwedge->otheredge->elabel = labels_g[i];
+
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			graph->removeEdge(spokes_g[i]);
+	};
+
+	PREPARE_REDUCTIONS2(m_g, spokes_g);
+	CHECK_REDUCTIONS(spokes_g);
+	H_DEBUG("Apply removeSpokeCenters");
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, spokes_g);
+	return true;
+}
+
+bool HopcroftImpl::removeStars()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<Vertex*> stars_g1;
+	std::vector<Vertex*> stars_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(stars_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isStar()) {
+				numbers.resize(numbers.size() + 1);
+				vertex->forEachEdge([&](EdgeEnd *e) {
+					EdgeEnd *f = e->cwedge;
+					numbers[ni].push_back(vertex->vlabel);
+					numbers[ni].push_back(e->elabel);
+					numbers[ni].push_back(f->elabel);
+					return true;
+				});
+				ni++;
+				stars_g.push_back(vertex);
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(stars_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			Vertex *star = stars_g[i];
+			star->forEachEdge([=](EdgeEnd *e) {
+				Vertex *v = e->isMultiEdge() ? nullptr : e->vertex;
+				graph->removeEdge(e, /* removeVertices */ false);
+				if (v)
+					graph->removeVertex(v);
+				return true;
+			});
+			star->vlabel = labels_g[i];
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, stars_g);
+	CHECK_REDUCTIONS(stars_g);
+	H_DEBUG("Apply removeStars");
+	CIRCLE_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, stars_g);
+	return true;
+}
+bool HopcroftImpl::removeDumbells()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<Vertex*> dumbells_g1;
+	std::vector<Vertex*> dumbells_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(dumbells_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->isDumbell()) {
+				numbers.resize(numbers.size() + 1);
+
+				numbers[ni].push_back(vertex->vlabel);
+				numbers[ni].push_back(vertex->vedge->elabel);
+				numbers[ni].push_back(vertex->vedge->otheredge->elabel);
+				numbers[ni].push_back(vertex->vlabel);
+				numbers[ni].push_back(vertex->vedge->otheredge->elabel);
+				numbers[ni].push_back(vertex->vedge->elabel);
+
+				ni++;
+				dumbells_g.push_back(vertex);
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(dumbells_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			Vertex *dumbell = dumbells_g[i];
+			if (dumbell->vedge) // eh ...
+				graph->removeEdge(dumbell->vedge, /* removeVertices */ false);
+			dumbell->vlabel = labels_g[i];
+			if (graph->vertices.size() > 1)
+				graph->removeVertex(dumbell);
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, dumbells_g);
+	CHECK_REDUCTIONS(dumbells_g);
+	H_DEBUG("Apply removeDumbells");
+	CIRCLE_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, dumbells_g);
+	return true;
+}
+
+bool HopcroftImpl::reductionB()
+{
+	// Bond associated reductions
+	if (replaceClumps())
+		return true;
+	if (replaceSkeins())
+		return true;
+	return false;
+}
+
+// B(i) Replacing clumps
+bool HopcroftImpl::replaceClumps()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<std::vector<EdgeEnd*>> clumps_g1;
+	std::vector<std::vector<EdgeEnd*>> clumps_g2;
+	std::vector<std::pair<Vertex*, Vertex*>> clumpvertices_g1;
+	std::vector<std::pair<Vertex*, Vertex*>> clumpvertices_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(clumps_g), ARG_REF(clumpvertices_g)) {
+		int ni = numbers.size();
+		for (EdgeEnd *edge : graph->edge_ends) {
+			Vertex *v = edge->vertex;
+			Vertex *w = edge->otheredge->vertex;
+
+			if (v->degree < 2 && w->degree < 2)
+				continue;
+
+			if (contains(clumpvertices_g, std::make_pair(v, w)) || contains(clumpvertices_g, std::make_pair(w, v)))
+				continue;
+
+			const std::vector<std::vector<EdgeEnd*>> &clumps = graph->findClumps(v, w);
+			if (clumps.empty())
+				continue;
+
+			clumps_g.insert(clumps_g.end(), clumps.begin(), clumps.end());
+			clumpvertices_g.push_back(std::make_pair(v, w));
+
+			for (const auto &clump : clumps) {
+				numbers.resize(numbers.size() + 2);
+
+				for (size_t i = 0; i < clump.size(); ++i)
+					numbers[ni].push_back(clump[i]->otheredge->elabel);
+
+				for (size_t i = clump.size(); i-- > 0; )
+					numbers[ni + 1].push_back(clump[i]->elabel);
+
+				ni += 2;
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(clumps_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); i += 2) {
+			const auto &clump = clumps_g[i / 2];
+			for (size_t j = 0; j < clump.size() - 1; ++j)
+				graph->removeEdge(clump[j]);
+
+			EdgeEnd *last = clump.back();
+			last->elabel = labels_g[i];
+			last->otheredge->elabel = labels_g[i + 1];
+		}
+	};
+
+	PREPARE_REDUCTIONS3(m_g, clumps_g, clumpvertices_g);
+	CHECK_REDUCTIONS(clumps_g);
+	H_DEBUG("Apply replaceClumps");
+	CIRCLE_PARTITIONING(numbers, labels); // FIXME: Vector in paper!
+	APPLY_REDUCTIONS2(m_g, clumps_g);
+	return true;
+}
+
+// B(ii) Replacing skeins
+bool HopcroftImpl::replaceSkeins()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers;
+	std::vector<std::pair<Vertex*, Vertex*>> skeinvertices_g1;
+	std::vector<std::pair<Vertex*, Vertex*>> skeinvertices_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(skeinvertices_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (!vertex->vedge)
+				continue;
+
+			Vertex *v = vertex;
+			Vertex *w = vertex->vedge->vertex;
+
+			if (v == w || (v->degree != 1 && w->degree != 1))
+				continue;
+
+			if (contains(skeinvertices_g, std::make_pair(v, w)) || contains(skeinvertices_g, std::make_pair(w, v)))
+				continue;
+
+			const std::vector<EdgeEnd*> &edges = v->neighborEdges();
+			if (edges.size() < 2)
+				continue;
+
+			skeinvertices_g.push_back(std::make_pair(v, w));
+			numbers.resize(numbers.size() + 2);
+
+			for (size_t i = 0; i < edges.size(); ++i) {
+				numbers[ni].push_back(v->vlabel);
+				numbers[ni].push_back(edges[i]->otheredge->elabel);
+				numbers[ni].push_back(edges[i]->elabel);
+				numbers[ni].push_back(w->vlabel);
+			}
+
+			for (size_t i = edges.size(); i-- > 0; ) {
+				numbers[ni + 1].push_back(w->vlabel);
+				numbers[ni + 1].push_back(edges[i]->elabel);
+				numbers[ni + 1].push_back(edges[i]->otheredge->elabel);
+				numbers[ni + 1].push_back(v->vlabel);
+			}
+
+			ni += 2;
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(skeinvertices_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); i += 2) {
+			const auto &skein = skeinvertices_g[i / 2];
+			skein.first->forEachEdge([=](EdgeEnd *e) {
+				graph->removeEdge(e, /* removeVertices */ false);
+				return true;
+			});
+			graph->removeVertex(skein.second);
+			skein.first->vlabel = std::min(labels_g[i], labels_g[i + 1]);
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, skeinvertices_g);
+	CHECK_REDUCTIONS(skeinvertices_g);
+	H_DEBUG("Apply replaceSkeins");
+	CIRCLE_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, skeinvertices_g);
+	return true;
+}
+
+bool HopcroftImpl::reductionC()
+{
+	// Four general classes of reductions and two special cases
+	if (removeIsolatedLowDegreeVertices())
+		return true;
+	if (removeNonregularGraphsWithNoIsolatedLowDegreeVertices())
+		return true;
+	if (removeLowDegreeVerticesWithSolitaryFaces())
+		return true;
+	if (removeTheLastGeneralCase())
+		return true;
+	return false;
+}
+
+// C (i) Removing isolated low degre vertices
+bool HopcroftImpl::removeIsolatedLowDegreeVertices()
+{
+	if (removeIsolatedLowDegree(2) || m_fail)
+		return true;
+	if (removeIsolatedLowDegree(3) || m_fail)
+		return true;
+	if (removeIsolatedLowDegree(4) || m_fail)
+		return true;
+	if (removeIsolatedLowDegree(5) || m_fail)
+		return true;
+	return false;
+}
+
+bool HopcroftImpl::removeIsolatedLowDegree(int degree)
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 3-tuple
+	std::vector<std::vector<std::pair<EdgeEnd*, EdgeEnd*>>> addededges_g1;
+	std::vector<std::vector<std::pair<EdgeEnd*, EdgeEnd*>>> addededges_g2;
+	std::vector<Vertex*> isolatedvertices_g1;
+	std::vector<Vertex*> isolatedvertices_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(addededges_g), ARG_REF(isolatedvertices_g)) {
+		int ni = numbers.size();
+		for (Vertex *vertex : graph->vertices) {
+			if (vertex->degree == degree && vertex->isDegreeIsolated()) {
+				std::vector<std::pair<EdgeEnd*, EdgeEnd*>> addededges;
+				const std::vector<EdgeEnd*> &edges = vertex->neighborEdges();
+				for (size_t i = 0; i < edges.size(); ++i) {
+					EdgeEnd *e1 = edges[i];
+					EdgeEnd *e2 = (i == edges.size() - 1) ? edges[0] : edges[i + 1];
+					if (e1->vertex == e2->vertex) // eh ...
+						continue;
+
+					numbers.resize(numbers.size() + 2);
+
+					numbers[ni] = {vertex->vlabel, e1->otheredge->elabel, e1->elabel};
+					numbers[ni + 1] = {e2->elabel, e2->otheredge->elabel, vertex->vlabel};
+					addededges.push_back(std::make_pair(e1, e2));
+
+					ni += 2;
+				}
+				addededges_g.push_back(addededges);
+				isolatedvertices_g.push_back(vertex);
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(addededges_g), ARG_CREF(isolatedvertices_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		int ni = 0;
+		for (size_t i = 0; i < isolatedvertices_g.size(); ++i) {
+			for (size_t j = 0; j < addededges_g[i].size(); ++j) {
+				EdgeEnd *e1 = addededges_g[i][j].first;
+				EdgeEnd *e2 = addededges_g[i][j].second;
+				EdgeEnd *enew = graph->addInnerEdge(e1, e2);
+				enew->elabel = labels_g[ni++];
+				enew->otheredge->elabel = labels_g[ni++];
+			}
+		}
+
+		for (Vertex *vertex : isolatedvertices_g) {
+			vertex->forEachEdge([=](EdgeEnd *e) {
+				graph->removeEdge(e, /* removeVertices */ false);
+				return true;
+			});
+			graph->removeVertex(vertex);
+		}
+	};
+
+	PREPARE_REDUCTIONS3(m_g, addededges_g, isolatedvertices_g);
+	CHECK_REDUCTIONS(addededges_g);
+	H_DEBUG("Apply removeIsolatedLowDegree " << degree);
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS3(m_g, addededges_g, isolatedvertices_g);
+	return true;
+}
+
+static bool isFourDegreeException(const std::vector<EdgeEnd*> &nghbrs)
+{
+	if (nghbrs.size() != 4)
+		return false;
+
+	// 4, x, 4, y
+	if (nghbrs[0]->vertex->degree == 4)
+		return nghbrs[1]->vertex->degree != 4
+			&& nghbrs[2]->vertex->degree == 4
+			&& nghbrs[3]->vertex->degree != 4;
+
+	// x, 4, y, 4
+	return nghbrs[1]->vertex->degree == 4
+		&& nghbrs[2]->vertex->degree != 4
+		&& nghbrs[3]->vertex->degree == 4;
+}
+
+// C (ii) Handling nonregular graphs with no isolated vertices of low degree
+bool HopcroftImpl::removeNonregularGraphsWithNoIsolatedLowDegreeVertices()
+{
+	if (removeIsolatedDegree(2) || m_fail)
+		return true;
+	if (removeIsolatedDegree(3) || m_fail)
+		return true;
+	if (removeIsolatedDegree(4) || m_fail)
+		return true;
+	if (removeIsolatedDegree(5) || m_fail)
+		return true;
+	return false;
+}
+
+bool HopcroftImpl::removeIsolatedDegree(int degree)
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 4-tuple
+	std::vector<EdgeEnd*> distinguishededges_g1;
+	std::vector<EdgeEnd*> distinguishededges_g2;
+
+	// We need to make sure to choose the same distinguished edge in both graphs
+	// This is just too complex and doesn't even work in case the G1 is CW-ordered
+	// and G2 is CCW-ordered (probably also in some other).
+	//
+	// There really is not much in paper as how to find the distinguished edge...
+#define H_REMOVE_ISOLATED_DEGREE_HEURISTIC 1
+
+#if H_REMOVE_ISOLATED_DEGREE_HEURISTIC
+	struct CacheItem {
+		int label;
+		int rlabel;
+		int degree;
+		std::vector<int> circle;
+	};
+	std::vector<CacheItem> distinguishedCache;
+
+	auto makeCircle = [](const std::vector<EdgeEnd*> &edges) {
+		std::vector<int> circle;
+		for (EdgeEnd *e : edges) {
+			circle.push_back(e->vertex->degree);
+			circle.push_back(e->otheredge->elabel);
+			circle.push_back(e->elabel);
+			circle.push_back(e->vertex->degree);
+		}
+		return circle;
+	};
+
+	auto degreeFromCache = [&](const std::vector<EdgeEnd*> &edges) {
+		const std::vector<int> &circle = makeCircle(edges);
+		std::vector<int> rcircle(circle.begin(), circle.end());
+		std::reverse(rcircle.begin(), rcircle.end());
+		for (const auto &c : distinguishedCache) {
+			if (circleIsomorphism(c.circle, circle))
+				return std::make_pair(c.label, c.degree);
+			if (circleIsomorphism(c.circle, rcircle))
+				return std::make_pair(c.rlabel, c.degree);
+		}
+		return std::make_pair(-1, -1);
+	};
+#endif
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(distinguishededges_g)) {
+		for (Vertex *v : graph->vertices) {
+			if (v->degree != degree)
+				continue;
+
+			const std::vector<EdgeEnd*> &nghbrs = v->neighborEdges(CounterClockWise);
+			if (isFourDegreeException(nghbrs))
+				continue;
+
+#if H_REMOVE_ISOLATED_DEGREE_HEURISTIC
+			EdgeEnd *distig = nullptr;
+			const auto &cached = degreeFromCache(nghbrs);
+			if (cached.first != -1) {
+				bool afterSame = false;
+				for (EdgeEnd *e : nghbrs) {
+					if (!distig && e->vertex->degree == cached.second && e->elabel == cached.first)
+						distig = e;
+					if (e->vertex->degree == degree)
+						afterSame = true;
+					else if (afterSame && e->vertex->degree == cached.second && e->elabel == cached.first) {
+						distig = e;
+						break;
+					}
+				}
+			}
+
+			if (!distig)  {
+				bool sameDegree = false;
+				for (EdgeEnd *e : nghbrs) {
+					if (e->vertex->degree == degree)
+						sameDegree = true;
+					else if (e->vertex->degree != degree)
+						distig = e;
+				}
+
+				if (!sameDegree || !distig)
+					distig = nullptr;
+				else
+					distinguishedCache.push_back({distig->elabel, distig->otheredge->elabel, distig->vertex->degree, makeCircle(nghbrs)});
+			}
+
+			if (!distig)
+				continue;
+#else
+
+			bool sameDegree = false;
+			EdgeEnd *distig = nullptr;
+			for (EdgeEnd *e : nghbrs) {
+				if (e->vertex->degree == degree) {
+					sameDegree = true;
+				} else if (e->vertex->degree != degree) {
+					distig = e;
+					break;
+				}
+			}
+
+			if (!sameDegree || !distig)
+				continue;
+#endif
+
+			EdgeEnd *inc = distig->otheredge->cwedge;
+			numbers.push_back({distig->elabel, distig->otheredge->elabel, v->vlabel, inc->otheredge->elabel});
+			distinguishededges_g.push_back(distig);
+			H_ASSERT(!contains(distinguishededges_g, distig->otheredge));
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(distinguishededges_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			EdgeEnd *distig = distinguishededges_g[i];
+			distig->otheredge->cwedge->otheredge->elabel = labels_g[i];
+			graph->removeDistinguishedEdgeAndVertex(distig);
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, distinguishededges_g);
+	CHECK_REDUCTIONS(distinguishededges_g);
+	H_DEBUG("Apply removeIsolatedDegree " << degree);
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, distinguishededges_g);
+	return true;
+}
+
+// C (iii) The first exception
+bool HopcroftImpl::removeFourDegreeException()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 3-tuple
+	std::vector<std::pair<Vertex*, EdgeEnd*>> circles_g1;
+	std::vector<std::pair<Vertex*, EdgeEnd*>> circles_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(circles_g)) {
+		for (Vertex *v : graph->vertices) {
+			if (v->degree != 4)
+				continue;
+
+			const std::vector<EdgeEnd*> &nghbrs = v->neighborEdges(CounterClockWise);
+			if (!isFourDegreeException(nghbrs))
+				continue;
+
+			int s = nghbrs[0]->vertex->degree == 4 ? 0 : 1;
+			circles_g.push_back(std::make_pair(v, nghbrs[s]));
+			H_ASSERT(!contains(circles_g, std::make_pair(nghbrs[s]->vertex, nghbrs[s]->otheredge)));
+
+			EdgeEnd *er = nghbrs[s];
+			EdgeEnd *es = nghbrs[s + 1];
+			EdgeEnd *et = nghbrs[s + 2];
+			EdgeEnd *eu = nghbrs[s + 3];
+
+			numbers.push_back({v->vlabel, es->otheredge->elabel, er->otheredge->elabel});
+			numbers.push_back({v->vlabel, eu->otheredge->elabel, et->otheredge->elabel});
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(circles_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); i += 2) {
+			const auto &circle = circles_g[i / 2];
+			graph->removeFourDegreeExceptionVertex(circle.first, circle.second);
+			circle.second->otheredge->elabel = labels_g[i];
+			circle.second->otheredge->cwedge->otheredge->cwedge->otheredge->elabel = labels_g[i + 1];
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, circles_g);
+	CHECK_REDUCTIONS(circles_g);
+	H_DEBUG("Apply fourDegreeException");
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, circles_g);
+	return true;
+}
+
+// C (iv) Handling regular graphs with vertices with solitary faces
+bool HopcroftImpl::removeLowDegreeVerticesWithSolitaryFaces()
+{
+	// d = 2  - special case (handled in testByDefinition)
+	if (removeVertexWithSolitaryFace(3) || m_fail)
+		return true;
+	if (removeVertexWithSolitaryFace(4) || m_fail)
+		return true;
+	if (removeVertexWithSolitaryFace(5) || m_fail)
+		return true;
+	// more than five face sides
+	if (removeVertexWithSolitaryFace(SPECIAL_X) || m_fail)
+		return true;
+	return false;
+}
+
+static int specialFaceLabel(Face *f)
+{
+	H_ASSERT(f->edges > 2);
+	return f->edges > 5 ? SPECIAL_X : f->edges;
+}
+
+// C (iv) Handling regular graphs with vertices with solitary faces
+bool HopcroftImpl::removeVertexWithSolitaryFace(int sides)
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 4-tuple
+	std::vector<std::vector<EdgeEnd*>> trees_g1;
+	std::vector<std::vector<EdgeEnd*>> trees_g2;
+	std::vector<std::vector<EdgeEnd*>> circles_g1;
+	std::vector<std::vector<EdgeEnd*>> circles_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(trees_g), ARG_REF(circles_g)) {
+		// Find distinguished edges + vertices
+		std::unordered_map<EdgeEnd*, EdgeEnd*> omap;
+		std::unordered_map<Vertex*, std::vector<EdgeEnd*>> rmap;
+		for (Vertex *v : graph->vertices) {
+			if (v->degree < 3)
+				continue;
+
+			// 4 possible face labels (3, 4, 5, x)
+			std::vector<int> counts(4, 0);
+			std::vector<EdgeEnd*> edges(4, nullptr);
+
+			v->forEachEdge([&](EdgeEnd *e) {
+				int idx = specialFaceLabel(e->ccwface);
+				idx = idx == SPECIAL_X ? 3 : idx - 3;
+				counts[idx]++;
+				edges[idx] = e->otheredge;
+				return true;
+			});
+
+			EdgeEnd *distig = nullptr;
+			for (size_t i = 0; i < 4; ++i)
+				if (counts[i] == 1 && specialFaceLabel(edges[i]->otheredge->ccwface) == sides)
+					distig = edges[i];
+
+			if (!distig)
+				continue;
+
+			rmap[v].push_back(distig);
+			rmap[distig->otheredge->vertex].push_back(distig->otheredge);
+
+			omap[distig] = distig;
+			omap[distig->otheredge] = distig;
+		}
+
+		if (rmap.empty())
+			return;
+
+		auto nextEdge = [](const std::vector<EdgeEnd*> &edges, Vertex *v, Vertex *w) {
+			EdgeEnd *next;
+			if (edges.front()->vertex == v || edges.front()->otheredge->vertex == v)
+				next = edges.back();
+			else
+				next = edges.front();
+			if (next->vertex == w)
+				next = next->otheredge;
+			return next;
+		};
+
+		// Build trees
+		bool found;
+		do {
+			found = false;
+			for (const auto &p : rmap) {
+				if (p.second.size() != 1)
+					continue;
+
+				std::vector<EdgeEnd*> tree;
+				if (p.second.front()->vertex == p.first)
+					tree.push_back(p.second.front()->otheredge);
+				else
+					tree.push_back(p.second.front());
+				rmap.erase(p.first);
+
+				Vertex *v = tree.back()->vertex;
+				while (rmap.find(v) != rmap.end()) {
+					const std::vector<EdgeEnd*> &edges = rmap[v];
+					if (edges.size() == 1) {
+						rmap.erase(v);
+						break;
+					}
+					tree.push_back(nextEdge(edges, tree.back()->otheredge->vertex, v));
+					rmap.erase(v);
+					v = tree.back()->vertex;
+				}
+				found = true;
+
+				// Use the original distinguished edge-end
+				for (size_t i = 0; i < tree.size(); ++i)
+					tree[i] = omap[tree[i]];
+
+				trees_g.push_back(tree);
+				break;
+			}
+		} while (found);
+
+		// Build tree vectors
+		for (std::vector<EdgeEnd*> &edges : trees_g) {
+			for (size_t i = 0; i < edges.size(); i++) {
+				EdgeEnd *distig = edges[i];
+				EdgeEnd *distig2 = distig->otheredge;
+				EdgeEnd *inc = distig2->otheredge->cwedge;
+				Vertex *v = distig2->otheredge->vertex;
+				numbers.push_back({distig2->elabel, distig2->otheredge->elabel, v->vlabel, inc->otheredge->elabel});
+			}
+		}
+
+		// Build circles
+		while (!rmap.empty()) {
+			for (const auto &p : rmap) {
+				std::vector<EdgeEnd*> circle;
+				circle.push_back(p.second.front()->otheredge);
+				rmap.erase(p.first);
+
+				Vertex *v = circle.back()->vertex;
+				while (v != circle.front()->otheredge->vertex) {
+					const std::vector<EdgeEnd*> &edges = rmap[v];
+					circle.push_back(nextEdge(edges, circle.back()->otheredge->vertex, v));
+					rmap.erase(v);
+					v = circle.back()->vertex;
+				}
+
+				circles_g.push_back(circle);
+				break;
+			}
+		};
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(trees_g), ARG_CREF(circles_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		// Trees
+		std::vector<EdgeEnd*> distigs;
+		for (const std::vector<EdgeEnd*> &edges : trees_g)
+			for (EdgeEnd *distig : edges)
+				distigs.push_back(distig);
+
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			distigs[i]->otheredge->cwedge->otheredge->elabel = labels_g[i];
+
+		for (size_t i = 0; i < labels_g.size(); ++i)
+			graph->removeDistinguishedEdgeAndVertex(distigs[i]->otheredge);
+
+		// Circles
+		int ulabel = newInteger();
+		for (const std::vector<EdgeEnd*> edges : circles_g) {
+			for (size_t i = 0; i < edges.size() - 1; ++i) {
+				EdgeEnd *distig = edges[i];
+				// Need to remove multi-edge with 2 last vertices of circle
+				if (distig->ccwedge->vertex == distig->otheredge->vertex) {
+					distig->vertex->vlabel = ulabel;
+					graph->removeEdge(distig->ccwedge);
+				} else if (distig->cwedge->vertex == distig->otheredge->vertex) {
+					distig->vertex->vlabel = ulabel;
+					graph->removeEdge(distig->cwedge);
+				}
+				graph->removeDistinguishedEdgeAndVertex(distig);
+			}
+		}
+	};
+
+	PREPARE_REDUCTIONS3(m_g, trees_g, circles_g);
+	CHECK_REDUCTIONS2(circles_g, trees_g);
+	H_DEBUG("Apply removeVertexWithSolitaryFace " << sides);
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS3(m_g, trees_g, circles_g);
+	return true;
+}
+
+static int anotherSpecialFaceLabel(Face *f)
+{
+	H_ASSERT(f->edges > 2);
+	return f->edges > 3 ? SPECIAL_X : f->edges;
+}
+
+static bool isAnotherFourDegreeException(const std::vector<EdgeEnd*> &nghbrs)
+{
+	if (nghbrs.size() != 4)
+		return false;
+
+	int l1 = anotherSpecialFaceLabel(nghbrs[0]->ccwface);
+	int l2 = anotherSpecialFaceLabel(nghbrs[1]->ccwface);
+	int l3 = anotherSpecialFaceLabel(nghbrs[2]->ccwface);
+	int l4 = anotherSpecialFaceLabel(nghbrs[3]->ccwface);
+
+	// 3, x, 3, x
+	if (l1 == 3)
+		return l2 == SPECIAL_X && l3 == 3 && l4 == SPECIAL_X;
+
+	// x, 3, x, 3
+	return l2 == 3 && l3 == SPECIAL_X && l4 == 3;
+}
+
+// C (v) The last general case
+bool HopcroftImpl::removeTheLastGeneralCase()
+{
+	std::vector<int> labels;
+	std::vector<std::vector<int>> numbers; // 4-tuple
+	std::vector<EdgeEnd*> distinguishededges_g1;
+	std::vector<EdgeEnd*> distinguishededges_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(distinguishededges_g)) {
+		for (Vertex *v : graph->vertices) {
+			if (v->degree < 4)
+				continue;
+
+			const std::vector<EdgeEnd*> &nghbrs = v->neighborEdges();
+			if (isAnotherFourDegreeException(nghbrs))
+				continue;
+
+			// 2 possible face labels (3, x)
+			std::vector<int> counts(2, 0);
+			std::vector<EdgeEnd*> edges(2, nullptr);
+
+			for (EdgeEnd *e : nghbrs) {
+				int idx = anotherSpecialFaceLabel(e->ccwface);
+				idx = idx == SPECIAL_X ? 1 : 0;
+				counts[idx]++;
+				edges[idx] = e->otheredge;
+			}
+
+			EdgeEnd *distig = nullptr;
+			if (counts[0] == 2 && counts[1] != 2)
+				distig = edges[0];
+			else if (counts[0] > 0 && counts[0] != 2 && counts[1] == 2)
+				distig = edges[1];
+
+			if (!distig)
+				continue;
+
+			EdgeEnd *inci = distig->otheredge->cwedge;
+			numbers.push_back({distig->elabel, distig->otheredge->elabel, v->vlabel, inci->otheredge->elabel});
+			distinguishededges_g.push_back(distig);
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(distinguishededges_g)) {
+		const std::vector<int> &labels_g = GET_LABELS(labels);
+		for (size_t i = 0; i < labels_g.size(); ++i) {
+			int label = labels_g[i];
+			EdgeEnd *distig = distinguishededges_g[i];
+			EdgeEnd *distigother = distig->otheredge;
+			int otherpos = indexOf(distinguishededges_g, distigother);
+			if (otherpos == -1) {
+				distig->otheredge->cwedge->otheredge->elabel = label;
+				graph->removeDistinguishedEdgeAndVertex(distig);
+			} else if (otherpos < int(i)) {
+				// Special case
+				int labelother = labels_g[otherpos];
+				distig->otheredge->cwedge->otheredge->elabel = label;
+				distigother->otheredge->cwedge->otheredge->elabel = labelother;
+				distig->vertex->vlabel = newInteger();
+				graph->removeDistinguishedEdgeAndVertex(distig);
+			}
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, distinguishededges_g);
+	CHECK_REDUCTIONS(distinguishededges_g);
+	H_DEBUG("Apply removeTheLastGeneralCase");
+	VECTOR_PARTITIONING(numbers, labels);
+	APPLY_REDUCTIONS2(m_g, distinguishededges_g);
+	return true;
+}
+
+// C (vi) Another special degree 4 case
+bool HopcroftImpl::removeAnotherSpecialDegreeFourCase()
+{
+	std::vector<EdgeEnd*> distinguishededges_g1;
+	std::vector<EdgeEnd*> distinguishededges_g2;
+
+	PREPARE_FUNCTION(Graph *graph, ARG_REF(distinguishededges_g)) {
+		for (Vertex *v : graph->vertices) {
+			if (v->degree != 4)
+				continue;
+
+			const std::vector<EdgeEnd*> &nghbrs = v->neighborEdges();
+			if (!isAnotherFourDegreeException(nghbrs))
+				continue;
+
+			for (EdgeEnd *e : nghbrs) {
+				if (e->ccwface->edges == 4) {
+					distinguishededges_g.push_back(e);
+				}
+			}
+		}
+	};
+
+	APPLY_FUNCTION(Graph *graph, ARG_CREF(distinguishededges_g)) {
+		for (EdgeEnd *e : distinguishededges_g) {
+			graph->addAndReconnectVertexIntoTriangularFace(e);
+			// Should the edges be relabelled?
+		}
+	};
+
+	PREPARE_REDUCTIONS2(m_g, distinguishededges_g);
+	CHECK_REDUCTIONS(distinguishededges_g);
+	H_DEBUG("Apply removeAnotherSpecialDegreeFourCase");
+	APPLY_REDUCTIONS2(m_g, distinguishededges_g);
+	return true;
+}
+
+bool HopcroftImpl::testByDefinition()
+{
+	H_ASSERT(!m_g1->vertices.empty());
+	H_ASSERT(!m_g2->vertices.empty());
+
+	if (isRegularDegreeTwoGraph())
+		return regularDegreeTwoGraphIsomorphism();
+
+	if (m_g1->vertices.size() == 1)
+		return m_g1->vertices.front()->vlabel == m_g2->vertices.front()->vlabel;
+
+	std::unordered_map<Vertex*, Vertex*> mapping;
+	return isomorphism_r(m_g1->vertices, m_g2->vertices, mapping);
+}
+
+bool HopcroftImpl::isRegularDegreeTwoGraph()
+{
+	for (Vertex *v : m_g1->vertices)
+		if (v->degree != 2)
+			return false;
+
+	return true;
+}
+
+// [C(iv)] Regular degree two graphs are n-gons which are handled by applying
+// the circle isomorphism directly
+//
+//  * Note that two circles (one clockwise, one counterclockwise) have to be
+//    formed for each n-gon in an analogous manner to skeins (its dual case).
+bool HopcroftImpl::regularDegreeTwoGraphIsomorphism()
+{
+	auto createCircle = [](Graph *graph, Direction d) {
+		std::vector<int> circle;
+		EdgeEnd *s = graph->edge_ends.front();
+		if (d == CounterClockWise)
+			s = s->otheredge;
+		EdgeEnd *e = s;
+		do {
+			circle.push_back(e->vertex->vlabel);
+			circle.push_back(e->elabel);
+			e = e->cwedge;
+		} while (e != s);
+		return circle;
+	};
+
+	const std::vector<int> &cw1 = createCircle(m_g1, ClockWise);
+	const std::vector<int> &ccw1 = createCircle(m_g1, CounterClockWise);
+
+	const std::vector<int> &cw2 = createCircle(m_g2, ClockWise);
+	const std::vector<int> &ccw2 = createCircle(m_g2, CounterClockWise);
+
+	if (circleIsomorphism(cw1, cw2) && circleIsomorphism(ccw1, ccw2))
+		return true;
+
+	if (circleIsomorphism(cw1, ccw2) && circleIsomorphism(ccw1, cw2))
+		return true;
+
+	return false;
+}
+
+// Recursive exhaustive matching isomorphism test
+bool HopcroftImpl::isomorphism_r(std::vector<Vertex*> vertices1, std::vector<Vertex*> vertices2, std::unordered_map<Vertex*, Vertex*> &mapping)
+{
+	if (vertices1.empty())
+		return isomorphism_mapping(mapping);
+
+	std::vector<Vertex*> vrtcs1 = vertices1;
+	Vertex *vertex = vertices1.back();
+	vrtcs1.pop_back();
+
+	for (size_t i = 0; i < vertices2.size(); ++i) {
+		std::vector<Vertex*> vrtcs2 = vertices2;
+		mapping.erase(vertex);
+		mapping.insert({vertex, vertices2[i]});
+		vrtcs2.erase(vrtcs2.begin() + i);
+		if (isomorphism_r(vrtcs1, vrtcs2, mapping))
+			return true;
+	}
+
+	return false;
+}
+
+bool HopcroftImpl::isomorphism_mapping(std::unordered_map<Vertex*, Vertex*> &mapping)
+{
+	auto createCircle = [](Vertex *v, Direction d) {
+		std::vector<int> circle;
+		v->forEachEdge([&](EdgeEnd *e) {
+			circle.push_back(e->elabel);
+			return true;
+		}, d);
+		return circle;
+	};
+
+	int cmp = 0; // 1 = cw-cw, 2 = cw-ccw
+
+	for (Vertex *v : m_g1->vertices) {
+		const std::vector<int> &cw1 = createCircle(v, ClockWise);
+		const std::vector<int> &ccw1 = createCircle(v, CounterClockWise);
+
+		const std::vector<int> &cw2 = createCircle(mapping[v], ClockWise);
+		const std::vector<int> &ccw2 = createCircle(mapping[v], CounterClockWise);
+
+		bool cwcw = circleIsomorphism(cw1, cw2) && circleIsomorphism(ccw1, ccw2);
+		bool cwccw = circleIsomorphism(cw1, ccw2) && circleIsomorphism(ccw1, cw2);
+
+		if (cmp == 1 && !cwcw)
+			return false;
+		else if (cmp == 2 && !cwccw)
+			return false;
+		else if (cmp == 0 && cwcw)
+			cmp = 1;
+		else if (cmp == 0 && cwccw)
+			cmp = 2;
+		else if (cmp == 0)
+			return false;
+	}
+
+	return true;
+}
+
+void HopcroftImpl::dumpGraphs()
+{
+	std::cout << "G1:";
+	m_g1->dump();
+
+	std::cout << "G2:";
+	m_g2->dump();
+}
+
+void HopcroftImpl::dumpLabels()
+{
+	std::cout << "G1:";
+	m_g1->dumpLabels();
+
+	std::cout << "G2:";
+	m_g2->dumpLabels();
+}
+
+} // namespace hopcroft
+
+} // namespace isomorphism
+
+} // namespace graph
diff --git a/alib2algo/src/graph/isomorphism/HopcroftImpl.h b/alib2algo/src/graph/isomorphism/HopcroftImpl.h
new file mode 100644
index 0000000000000000000000000000000000000000..ae7ac158c75652000202a30b0943c5bbcb7c51a0
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftImpl.h
@@ -0,0 +1,91 @@
+#ifndef GRAPH_HOPCROFT_IMPL_H_
+#define GRAPH_HOPCROFT_IMPL_H_
+
+#include <list>
+#include <vector>
+#include <unordered_map>
+
+#include "HopcroftGraph.h"
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+namespace hopcroft
+{
+
+class HopcroftImpl
+{
+public:
+	explicit HopcroftImpl(Graph *g1, Graph *g2);
+
+	bool testIsomorphism();
+
+//private:
+
+	// Subalgorithms
+	int newInteger();
+	void saveInteger();
+	void restoreInteger();
+	std::vector<int> vectorPartitioning(const std::vector<std::vector<int>> &tuples);
+	bool circleIsomorphism(std::vector<int> a, std::vector<int> b);
+	std::vector<int> circlePartitioning(const std::vector<std::vector<int>> &circles);
+
+	// Reductions
+	bool doReductions();
+
+	// A
+	bool reductionA();
+	bool removeLoops();
+	bool removeCorollas();
+	bool removeKnots();
+	bool removeDegreeOneVertices();
+	bool removeSpokeCenters();
+	bool removeStars();
+	bool removeDumbells();
+
+	// B
+	bool reductionB();
+	bool replaceClumps();
+	bool replaceSkeins();
+
+	// C
+	bool reductionC();
+	bool removeIsolatedLowDegreeVertices();
+	bool removeIsolatedLowDegree(int degree);
+	bool removeNonregularGraphsWithNoIsolatedLowDegreeVertices();
+	bool removeIsolatedDegree(int degree);
+	bool removeFourDegreeException();
+	bool removeLowDegreeVerticesWithSolitaryFaces();
+	bool removeVertexWithSolitaryFace(int sides);
+	bool removeTheLastGeneralCase();
+	bool removeAnotherSpecialDegreeFourCase();
+
+	// Final
+	bool testByDefinition();
+
+	bool isRegularDegreeTwoGraph();
+	bool regularDegreeTwoGraphIsomorphism();
+
+	bool isomorphism_r(std::vector<Vertex*> vertices1, std::vector<Vertex*> vertices2, std::unordered_map<Vertex*, Vertex*> &mapping);
+	bool isomorphism_mapping(std::unordered_map<Vertex*, Vertex*> &mapping);
+
+	void dumpGraphs();
+	void dumpLabels();
+
+	Graph *m_g1;
+	Graph *m_g2;
+	bool m_fail;
+	int m_integer;
+	int m_savedInteger;
+};
+
+} // namespace hopcroft
+
+} // namespace isomorphism
+
+} // namespace graph
+
+#endif // GRAPH_HOPCROFT_IMPL_H_
diff --git a/alib2algo/src/graph/isomorphism/HopcroftWong.cpp b/alib2algo/src/graph/isomorphism/HopcroftWong.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..22bd63460b8040e5149a7815a7381942a8311c76
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftWong.cpp
@@ -0,0 +1,28 @@
+#include "HopcroftWong.h"
+#include "HopcroftGraph.h"
+#include "HopcroftImpl.h"
+
+#include <exception/AlibException.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+bool HopcroftWong::hopcroftwong(const Graph &first, const Graph &second)
+{
+	return getInstance().dispatch(first.getData(), second.getData());
+}
+
+bool HopcroftWong::hopcroftwong(const UndirectedGraph &first, const UndirectedGraph &second)
+{
+	hopcroft::Graph g1(first);
+	hopcroft::Graph g2(second);
+	hopcroft::HopcroftImpl impl(&g1, &g2);
+	return impl.testIsomorphism();
+}
+
+} // namespace isomorphism
+
+} // namespace graph
diff --git a/alib2algo/src/graph/isomorphism/HopcroftWong.h b/alib2algo/src/graph/isomorphism/HopcroftWong.h
new file mode 100644
index 0000000000000000000000000000000000000000..85ac2a546cbbef948f8f597c4af2a689b2927c83
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/HopcroftWong.h
@@ -0,0 +1,34 @@
+#ifndef GRAPH_HOPCROFT_WONG_H_
+#define GRAPH_HOPCROFT_WONG_H_
+
+#include <common/multipleDispatch.hpp>
+
+#include <graph/Graph.h>
+#include <graph/directed/DirectedGraph.h>
+#include <graph/undirected/UndirectedGraph.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+// Decides isomorphism of triconnected planar graphs
+class HopcroftWong : public std::DoubleDispatch<bool, graph::GraphBase, graph::GraphBase>
+{
+public:
+	static bool hopcroftwong(const Graph &first, const Graph &second);
+
+	static bool hopcroftwong(const UndirectedGraph &first, const UndirectedGraph &second);
+
+	static HopcroftWong &getInstance() {
+		static HopcroftWong res;
+		return res;
+	}
+};
+
+} // namespace isomorphism
+
+} // namespace graph
+
+#endif // GRAPH_HOPCROFT_WONG_H_
diff --git a/alib2algo/src/graph/isomorphism/Isomorphism.cpp b/alib2algo/src/graph/isomorphism/Isomorphism.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5009e3f1d5798e3035760a93da1617b62d9c8535
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/Isomorphism.cpp
@@ -0,0 +1,124 @@
+#include "Isomorphism.h"
+
+#include <exception/AlibException.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+static inline std::set<UndirectedEdge> stripNames(const std::set<UndirectedEdge> &s)
+{
+	std::set<UndirectedEdge> out;
+	for (const UndirectedEdge &v : s)
+		out.insert(UndirectedEdge(v.getFirstNode(), v.getSecondNode()));
+	return out;
+}
+
+static inline std::set<DirectedEdge> stripNames(const std::set<DirectedEdge> &s)
+{
+	std::set<DirectedEdge> out;
+	for (const DirectedEdge &v : s)
+		out.insert(DirectedEdge(v.getFromNode(), v.getToNode()));
+	return out;
+}
+
+static bool is_equal(const std::set<UndirectedEdge> &first, const std::set<UndirectedEdge> &second, const std::unordered_map<Node, Node> &mapping)
+{
+	if (first.size() != second.size())
+		return false;
+
+	std::set<UndirectedEdge> frst = stripNames(first);
+	std::set<UndirectedEdge> scnd = stripNames(second);
+
+	for (const UndirectedEdge &e : frst) {
+		UndirectedEdge other1(mapping.at(e.getFirstNode()), mapping.at(e.getSecondNode()), e.getName());
+		UndirectedEdge other2(mapping.at(e.getSecondNode()), mapping.at(e.getFirstNode()), e.getName());
+		if (scnd.find(other1) == scnd.end() && scnd.find(other2) == scnd.end())
+			return false;
+	}
+
+	return true;
+}
+
+static bool is_equal(const std::set<DirectedEdge> &first, const std::set<DirectedEdge> &second, const std::unordered_map<Node, Node> &mapping)
+{
+	if (first.size() != second.size())
+		return false;
+
+	std::set<DirectedEdge> frst = stripNames(first);
+	std::set<DirectedEdge> scnd = stripNames(second);
+
+	for (const DirectedEdge &e : frst) {
+		DirectedEdge other(mapping.at(e.getFromNode()), mapping.at(e.getToNode()), e.getName());
+		if (scnd.find(other) == scnd.end())
+			return false;
+	}
+
+	return true;
+}
+
+template<typename T>
+static bool impl_r(const T &first, const T &second, const std::vector<Node> &nodes1, const std::vector<Node> &nodes2, std::unordered_map<Node, Node> &mapping)
+{
+	if (nodes1.empty()) {
+		for (const Node &n : first.getNodes()) {
+			if (!is_equal(first.neighborEdges(n), second.neighborEdges(mapping.at(n)), mapping)) {
+				return false;
+			}
+		}
+		return true;
+	}
+
+	std::vector<Node> nds1 = nodes1;
+	Node node = nodes1.back();
+	nds1.pop_back();
+
+	for (size_t i = 0; i < nodes2.size(); ++i) {
+		std::vector<Node> nds2 = nodes2;
+		mapping.erase(node);
+		mapping.insert({node, nodes2[i]});
+		nds2.erase(nds2.begin() + i);
+		if (impl_r(first, second, nds1, nds2, mapping))
+			return true;
+	}
+
+	return false;
+}
+
+template<typename T>
+static bool impl(const T &first, const T &second)
+{
+	const std::set<Node> &nodes1 = first.getNodes();
+	const std::set<Node> &nodes2 = second.getNodes();
+
+	if (nodes1.size() != nodes2.size())
+		return false;
+
+	std::vector<Node> n1;
+	std::vector<Node> n2;
+	std::unordered_map<Node, Node> mapping;
+	std::copy(nodes1.begin(), nodes1.end(), std::back_inserter(n1));
+	std::copy(nodes2.begin(), nodes2.end(), std::back_inserter(n2));
+	return impl_r(first, second, n1, n2, mapping);
+}
+
+bool Isomorphism::isomorphism(const Graph &first, const Graph &second)
+{
+	return getInstance().dispatch(first.getData(), second.getData());
+}
+
+bool Isomorphism::isomorphism(const DirectedGraph &first, const DirectedGraph &second)
+{
+	return impl(first, second);
+}
+
+bool Isomorphism::isomorphism(const UndirectedGraph &first, const UndirectedGraph &second)
+{
+	return impl(first, second);
+}
+
+} // namespace isomorphism
+
+} // namespace graph
diff --git a/alib2algo/src/graph/isomorphism/Isomorphism.h b/alib2algo/src/graph/isomorphism/Isomorphism.h
new file mode 100644
index 0000000000000000000000000000000000000000..c99499d4c53f3ee0e85426ed38186478923eee50
--- /dev/null
+++ b/alib2algo/src/graph/isomorphism/Isomorphism.h
@@ -0,0 +1,34 @@
+#ifndef GRAPH_ISOMORPHISM_H_
+#define GRAPH_ISOMORPHISM_H_
+
+#include <common/multipleDispatch.hpp>
+
+#include <graph/Graph.h>
+#include <graph/directed/DirectedGraph.h>
+#include <graph/undirected/UndirectedGraph.h>
+
+namespace graph
+{
+
+namespace isomorphism
+{
+
+class Isomorphism : public std::DoubleDispatch<bool, graph::GraphBase, graph::GraphBase>
+{
+public:
+	static bool isomorphism(const Graph &first, const Graph &second);
+
+	static bool isomorphism(const DirectedGraph &first, const DirectedGraph &second);
+	static bool isomorphism(const UndirectedGraph &first, const UndirectedGraph &second);
+
+	static Isomorphism &getInstance() {
+		static Isomorphism res;
+		return res;
+	}
+};
+
+} // namespace isomorphism
+
+} // namespace graph
+
+#endif // GRAPH_ISOMORPHISM_H_
diff --git a/alib2algo/test-src/graph/isomorphism/HopcroftTest.cpp b/alib2algo/test-src/graph/isomorphism/HopcroftTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c1e0291efd890a24931f1b6db69a1eff5fc16c24
--- /dev/null
+++ b/alib2algo/test-src/graph/isomorphism/HopcroftTest.cpp
@@ -0,0 +1,804 @@
+#include "HopcroftTest.h"
+
+#include <algorithm>
+
+#include "graph/isomorphism/HopcroftImpl.h"
+#include "graph/embedding/HopcroftTarjan.h"
+#include "graph/generate/RandomGraphFactory.h"
+
+CPPUNIT_TEST_SUITE_REGISTRATION(GraphHopcroftTest);
+CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(GraphHopcroftTest, "graph");
+
+static graph::UndirectedGraph getGraph1()
+{
+	graph::Node n1("n1");
+	graph::Node n2("n2");
+	graph::Node n3("n3");
+	graph::Node n4("n4");
+	graph::Node n5("n5");
+
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+	ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+
+	return ug1;
+}
+
+static graph::UndirectedGraph getGraph1_I()
+{
+	graph::Node n1("ni1i");
+	graph::Node n2("ni2i");
+	graph::Node n3("ni3i");
+	graph::Node n4("ni4i");
+	graph::Node n5("ni5i");
+
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n2, n1));
+	ug1.addEdge(graph::UndirectedEdge(n3, n2));
+	ug1.addEdge(graph::UndirectedEdge(n4, n1));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n3, n5));
+	ug1.addEdge(graph::UndirectedEdge(n3, n1));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+	ug1.addEdge(graph::UndirectedEdge(n4, n3));
+	ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+
+	return ug1;
+}
+
+static graph::UndirectedGraph getGraph2()
+{
+	graph::Node n1("ni1");
+	graph::Node n2("ni2");
+	graph::Node n3("ni3");
+	graph::Node n4("ni4");
+	graph::Node n5("ni5");
+
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge1"));
+	ug1.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge2"));
+	ug1.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge3"));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4, "multi-edge4"));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3, "multi-edge5"));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4, "multi-edge6"));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+
+	std::unordered_map<graph::Node, std::vector<graph::Node>> embedding;
+	embedding[n1] = { n4, n2, n2, n2, n2, n3, n4 };
+	embedding[n2] = { n1, n1, n5, n3, n3, n1, n1 };
+	embedding[n3] = { n4, n1, n2, n2, n5 };
+	embedding[n4] = { n1, n1, n3, n5, n5 };
+	embedding[n5] = { n3, n2, n4, n4 };
+	ug1.setEmbedding(embedding);
+
+	return ug1;
+}
+
+static graph::UndirectedGraph getGraph2_I()
+{
+	graph::Node n1("n1i");
+	graph::Node n2("n2i");
+	graph::Node n3("n3i");
+	graph::Node n4("n4i");
+	graph::Node n5("n5i");
+
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n2, n1, "multi-edge3"));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n4, n1, "multi-edge4"));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n2, n1, "multi-edge2"));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3, "multi-edge5"));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n4, n5));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4, "multi-edge6"));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge1"));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+
+	std::unordered_map<graph::Node, std::vector<graph::Node>> embedding;
+	embedding[n1] = { n3, n4, n4, n2, n2, n2, n2 };
+	embedding[n2] = { n5, n3, n3, n1, n1, n1, n1 };
+	embedding[n3] = { n5, n4, n1, n2, n2 };
+	embedding[n4] = { n3, n5, n5, n1, n1 };
+	embedding[n5] = { n3, n2, n4, n4 };
+	ug1.setEmbedding(embedding);
+
+	return ug1;
+}
+
+void GraphHopcroftTest::testCircleIsomorphism()
+{
+	{
+		std::vector<int> a = { 1, 2, 3, 4, 5 };
+		std::vector<int> b = { 1, 2, 3, 4, 5 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(true, impl.circleIsomorphism(a, b));
+	}
+
+	{
+		std::vector<int> a = { 1, 2, 3, 4, 5 };
+		std::vector<int> b = { 2, 3, 4, 5, 1 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(true, impl.circleIsomorphism(a, b));
+	}
+
+	{
+		std::vector<int> a = { 1, 2, 3, 4, 5 };
+		std::vector<int> b = { 5, 1, 2, 3, 4 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(true, impl.circleIsomorphism(a, b));
+	}
+
+	{
+		std::vector<int> a = { 1, 2, 3, 4, 5 };
+		std::vector<int> b = { 3, 4, 5, 1, 2 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(true, impl.circleIsomorphism(a, b));
+	}
+
+	{
+		std::vector<int> a = { 1, 2, 4, 4, 5 };
+		std::vector<int> b = { 2, 3, 4, 5, 1 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(false, impl.circleIsomorphism(a, b));
+	}
+
+	{
+		std::vector<int> a = { 1, 2, 4, 5 };
+		std::vector<int> b = { 2, 3, 4, 5, 1 };
+		graph::isomorphism::hopcroft::HopcroftImpl impl(nullptr, nullptr);
+		CPPUNIT_ASSERT_EQUAL(false, impl.circleIsomorphism(a, b));
+	}
+}
+
+static const int MAX_STEPS = 1000;
+
+#define RET_ERROR(info, expected, got) \
+	{ \
+		std::cout << info << " [ Expected:" << expected << " Got:" << got << " ]" << std::endl; \
+		return false; \
+	}
+
+static graph::isomorphism::hopcroft::Vertex *findVertex(graph::isomorphism::hopcroft::Graph *g, const std::string &l)
+{
+	for (graph::isomorphism::hopcroft::Vertex *v : g->vertices)
+		if (v->out == l)
+			return v;
+
+	CPPUNIT_ASSERT(false && "vertex not found");
+	return nullptr;
+}
+
+static graph::isomorphism::hopcroft::EdgeEnd *findEdgeEnd(graph::isomorphism::hopcroft::Graph *g, const std::string &l)
+{
+	for (graph::isomorphism::hopcroft::EdgeEnd *e : g->edge_ends)
+		if (l == e->otheredge->vertex->out + "->" + e->vertex->out)
+			return e;
+
+	CPPUNIT_ASSERT(false && "edge not found");
+	return nullptr;
+}
+
+static bool checkEdge(graph::isomorphism::hopcroft::EdgeEnd *e, graph::isomorphism::hopcroft::Graph *g)
+{
+	if (e->isSpoke()) {
+
+	}
+	if (e->isLoopEdge()) {
+		if (!e->otheredge->isLoopEdge())
+			RET_ERROR("e->otheredge not loop edge", true, false);
+	}
+	if (e->isMultiEdge()) {
+		if (!e->otheredge->isMultiEdge())
+			RET_ERROR("e->otheredge not multi edge", true, false);
+	}
+
+	if (std::find(g->edge_ends.begin(), g->edge_ends.end(), e->cwedge) == g->edge_ends.end())
+		RET_ERROR("cwedge not in graph", false, false);
+
+	if (std::find(g->edge_ends.begin(), g->edge_ends.end(), e->ccwedge) == g->edge_ends.end())
+		RET_ERROR("ccwedge not in graph", false, false);
+
+	if (std::find(g->edge_ends.begin(), g->edge_ends.end(), e->otheredge) == g->edge_ends.end())
+		RET_ERROR("otheredge not in graph", false, false);
+
+	if (e->ccwface && std::find(g->faces.begin(), g->faces.end(), e->ccwface) == g->faces.end())
+		RET_ERROR("ccwface not in graph", false, false);
+
+	if (std::find(g->vertices.begin(), g->vertices.end(), e->vertex) == g->vertices.end())
+		RET_ERROR("vertex not in graph", false, false);
+
+	return true;
+}
+
+static bool checkIncidentEdges(graph::isomorphism::hopcroft::Vertex *v)
+{
+	int i = 0;
+	std::set<graph::isomorphism::hopcroft::Vertex*> neighbors;
+	v->forEachEdge([&i, &neighbors](graph::isomorphism::hopcroft::EdgeEnd *e) {
+		neighbors.insert(e->vertex);
+		return ++i <= MAX_STEPS;
+	});
+
+	if (i > MAX_STEPS)
+		RET_ERROR("Over max_steps clockwise", false, false);
+
+	if ((int)neighbors.size() != v->degree)
+		RET_ERROR("Neighbors not equal to degree", neighbors.size(), v->degree);
+
+	i = 0;
+	v->forEachEdge([&i](graph::isomorphism::hopcroft::EdgeEnd *) {
+		return ++i <= MAX_STEPS;
+	});
+
+	if (i > MAX_STEPS)
+		RET_ERROR("Over max_steps counter-clockwise", false, false);
+
+	return true;
+}
+
+static bool checkVEdge(graph::isomorphism::hopcroft::Vertex *v, graph::isomorphism::hopcroft::Graph *g)
+{
+	if (v->degree == 0 && !v->vedge)
+		return true;
+
+	if (std::find(g->edge_ends.begin(), g->edge_ends.end(), v->vedge) == g->edge_ends.end())
+		RET_ERROR("vedge not in graph", false, false);
+
+	int i = 0;
+	graph::isomorphism::hopcroft::EdgeEnd *e = v->vedge;
+	while (i == 0 || e != v->vedge) {
+		if (++i > MAX_STEPS)
+			RET_ERROR("Over max_steps", false, false);
+		e = e->ccwedge;
+	}
+
+	return true;
+}
+
+static bool checkFace(graph::isomorphism::hopcroft::Face *f, graph::isomorphism::hopcroft::Graph *g)
+{
+	if (g->faces.size() == 1 && g->edge_ends.empty())
+		return true;
+
+	int i = 0;
+	graph::isomorphism::hopcroft::EdgeEnd *e = f->fedge;
+
+	while (i == 0 || e != f->fedge) {
+		if (++i > MAX_STEPS)
+			RET_ERROR("Over max_steps", false, false);
+		if (e->ccwface != f)
+			RET_ERROR("e->ccwface not f", f, e->ccwface);
+		e = e->ccwedge;
+	}
+
+	if (i != f->edges)
+		RET_ERROR("incorrect f->edges", i, f->edges);
+
+	if (std::find(g->edge_ends.begin(), g->edge_ends.end(), f->fedge) == g->edge_ends.end())
+		RET_ERROR("fedge not in graph", false, false);
+
+	return true;
+}
+
+static void checkVisitAllEdges(graph::isomorphism::hopcroft::EdgeEnd *e, graph::isomorphism::hopcroft::Graph *g, std::set<graph::isomorphism::hopcroft::EdgeEnd*> &edges, std::set<graph::isomorphism::hopcroft::Vertex*> &vertices)
+{
+	if (edges.find(e) != edges.end())
+		return;
+
+	edges.insert(e);
+	vertices.insert(e->vertex);
+
+	checkEdge(e, g);
+	checkVisitAllEdges(e->otheredge, g, edges, vertices);
+	checkVisitAllEdges(e->cwedge, g, edges, vertices);
+	checkVisitAllEdges(e->ccwedge, g, edges, vertices);
+}
+
+static bool checkGraphConsistency(graph::isomorphism::hopcroft::Graph *g)
+{
+	for (graph::isomorphism::hopcroft::Vertex *v : g->vertices) {
+		if (!checkVEdge(v, g)) {
+			std::cout << "FAIL: checkVEdge: " << v << std::endl;
+			return false;
+		}
+		if (!checkIncidentEdges(v)) {
+			std::cout << "FAIL: checkIncidentEdges: " << v << std::endl;
+			return false;
+		}
+	}
+
+	for (graph::isomorphism::hopcroft::EdgeEnd *e : g->edge_ends) {
+		if (!checkEdge(e, g)) {
+			std::cout << "FAIL: checkEdge: " << e << std::endl;
+			return false;
+		}
+	}
+
+	for (graph::isomorphism::hopcroft::Face *f : g->faces) {
+		if (!checkFace(f, g)) {
+			std::cout << "FAIL: checkFace: " << f << std::endl;
+			return false;
+		}
+	}
+
+	if (!g->edge_ends.empty()) {
+		std::set<graph::isomorphism::hopcroft::EdgeEnd*> edges;
+		std::set<graph::isomorphism::hopcroft::Vertex*> vertices;
+		checkVisitAllEdges(g->edge_ends.front(), g, edges, vertices);
+		if (edges.size() != g->edge_ends.size())
+			return false;
+	}
+	return true;
+}
+
+void GraphHopcroftTest::testEmbedding()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::isomorphism::hopcroft::Graph g1(ug1);
+	CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+	graph::UndirectedGraph ug2 = getGraph2();
+	graph::isomorphism::hopcroft::Graph g2(ug2);
+	CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+}
+
+void GraphHopcroftTest::testEdges()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::isomorphism::hopcroft::Graph g1(ug1);
+
+	CPPUNIT_ASSERT(checkIncidentEdges(findVertex(&g1, "n1")));
+	CPPUNIT_ASSERT(checkIncidentEdges(findVertex(&g1, "n2")));
+	CPPUNIT_ASSERT(checkIncidentEdges(findVertex(&g1, "n3")));
+	CPPUNIT_ASSERT(checkIncidentEdges(findVertex(&g1, "n4")));
+	CPPUNIT_ASSERT(checkIncidentEdges(findVertex(&g1, "n5")));
+}
+
+void GraphHopcroftTest::testRemoveEdge()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::UndirectedGraph ug2 = getGraph2();
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		CPPUNIT_ASSERT_EQUAL(5, (int)g1.faces.size());
+		CPPUNIT_ASSERT_EQUAL(5, (int)g1.vertices.size());
+		CPPUNIT_ASSERT_EQUAL(16, (int)g1.edge_ends.size());
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		while (!g1.edge_ends.empty()) {
+			g1.removeEdge(g1.edge_ends.front());
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+		}
+
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.faces.size());
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.vertices.size());
+		for (graph::isomorphism::hopcroft::Vertex *v : g1.vertices)
+			CPPUNIT_ASSERT(v->vedge == nullptr);
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		while (!g1.edge_ends.empty()) {
+			g1.removeEdge(g1.edge_ends.back());
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+		}
+
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.faces.size());
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.vertices.size());
+		for (graph::isomorphism::hopcroft::Vertex *v : g1.vertices)
+			CPPUNIT_ASSERT(v->vedge == nullptr);
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		while (!g1.vertices.empty()) {
+			graph::isomorphism::hopcroft::Vertex *v = g1.vertices.front();
+			v->forEachEdge([&g1](graph::isomorphism::hopcroft::EdgeEnd *e) {
+				g1.removeEdge(e, /* removeVertices */ false);
+				CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+				return true;
+			});
+
+			CPPUNIT_ASSERT_EQUAL(0, v->degree);
+			g1.removeVertex(v);
+		}
+
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.faces.size());
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.edge_ends.size());
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		while (!g1.vertices.empty()) {
+			graph::isomorphism::hopcroft::Vertex *v = g1.vertices.back();
+			v->forEachEdge([&g1](graph::isomorphism::hopcroft::EdgeEnd *e) {
+				g1.removeEdge(e, /* removeVertices */ false);
+				CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+				return true;
+			});
+
+			CPPUNIT_ASSERT_EQUAL(0, v->degree);
+			g1.removeVertex(v);
+		}
+
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.faces.size());
+		CPPUNIT_ASSERT_EQUAL(0, (int)g1.edge_ends.size());
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g2(ug2);
+
+		while (!g2.vertices.empty()) {
+			graph::isomorphism::hopcroft::Vertex *v = g2.vertices.back();
+			v->forEachEdge([&g2](graph::isomorphism::hopcroft::EdgeEnd *e) {
+				g2.removeEdge(e, /* removeVertices */ false);
+				CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+				return true;
+			});
+
+			CPPUNIT_ASSERT_EQUAL(0, v->degree);
+			g2.removeVertex(v);
+		}
+
+		CPPUNIT_ASSERT_EQUAL(0, (int)g2.faces.size());
+		CPPUNIT_ASSERT_EQUAL(0, (int)g2.edge_ends.size());
+	}
+}
+
+void GraphHopcroftTest::testAddInnerEdge()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		for (graph::isomorphism::hopcroft::Vertex *v : g1.vertices) {
+			std::vector<graph::isomorphism::hopcroft::EdgeEnd*> edges;
+			v->forEachEdge([&edges](graph::isomorphism::hopcroft::EdgeEnd *e) {
+				edges.push_back(e);
+				return true;
+			});
+
+			for (size_t i = 0; i < edges.size(); ++i) {
+				auto *e1 = edges[i];
+				auto *e2 = (i == edges.size() - 1) ? edges[0] : edges[i + 1];
+				if (e1->vertex == e2->vertex)
+					continue;
+
+				g1.addInnerEdge(e1, e2);
+				CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			}
+		}
+
+		while (!g1.edge_ends.empty()) {
+			g1.removeEdge(g1.edge_ends.back());
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+		}
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		g1.removeEdge(findEdgeEnd(&g1, "n2->n5"));
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n2->n1"));
+
+		g1.addInnerEdge(findEdgeEnd(&g1, "n3->n2"), findEdgeEnd(&g1, "n3->n5"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n2"));
+		g1.removeEdge(findEdgeEnd(&g1, "n3->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n4->n1"));
+
+		g1.addInnerEdge(findEdgeEnd(&g1, "n3->n2"), findEdgeEnd(&g1, "n3->n5"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+	}
+}
+
+void GraphHopcroftTest::testRemoveDistinguishedEdge()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		g1.removeDistinguishedEdgeAndVertex(findEdgeEnd(&g1, "n3->n4"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		g1.removeEdge(findEdgeEnd(&g1, "n1->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n2->n5"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		g1.removeDistinguishedEdgeAndVertex(findEdgeEnd(&g1, "n2->n4"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n4->n1"));
+		g1.removeEdge(findEdgeEnd(&g1, "n4->n1"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+	}
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n2"));
+		g1.removeEdge(findEdgeEnd(&g1, "n2->n1"));
+		g1.removeEdge(findEdgeEnd(&g1, "n2->n3"));
+		g1.removeEdge(findEdgeEnd(&g1, "n5->n4"));
+		g1.removeEdge(findEdgeEnd(&g1, "n1->n4"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+
+		g1.removeDistinguishedEdgeAndVertex(findEdgeEnd(&g1, "n3->n4"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+	}
+}
+
+void GraphHopcroftTest::testRemoveFourDegreeExceptionEdge()
+{
+	graph::Node n1("n1");
+	graph::Node n2("n2");
+	graph::Node n3("n3");
+	graph::Node n4("n4");
+	graph::Node n5("n5");
+
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n2, n4));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+	ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		g1.removeFourDegreeExceptionVertex(findVertex(&g1, "n3"), findEdgeEnd(&g1, "n3->n2"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+	}
+}
+
+void GraphHopcroftTest::testAddAndReconnectVertexIntoTriangularFace()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+
+	{
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+
+		g1.addAndReconnectVertexIntoTriangularFace(findEdgeEnd(&g1, "n3->n2"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+		g1.addAndReconnectVertexIntoTriangularFace(findEdgeEnd(&g1, "n3->n4"));
+		CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+		CPPUNIT_ASSERT_EQUAL(7, (int)g1.vertices.size());
+	}
+}
+
+void GraphHopcroftTest::testIsomorphismByDefinition()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::UndirectedGraph ug1_i = getGraph1_I();
+	graph::UndirectedGraph ug2 = getGraph2();
+	graph::UndirectedGraph ug2_i = getGraph2_I();
+
+	{
+		// Same simple graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(true, impl.testByDefinition());
+	}
+
+	{
+		// Same graphs with multi-edges and loops
+		graph::isomorphism::hopcroft::Graph g1(ug2);
+		graph::isomorphism::hopcroft::Graph g2(ug2);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(true, impl.testByDefinition());
+	}
+
+	{
+		// Isomorphic simple graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1_i);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(true, impl.testByDefinition());
+	}
+
+	{
+		// Isomorphic graphs with multi-edges and loops
+		graph::isomorphism::hopcroft::Graph g1(ug2);
+		graph::isomorphism::hopcroft::Graph g2(ug2_i);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(true, impl.testByDefinition());
+	}
+
+	{
+		// Different graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug2);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(false, impl.testByDefinition());
+	}
+
+	{
+		// Different graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug2_i);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		CPPUNIT_ASSERT_EQUAL(false, impl.testByDefinition());
+	}
+}
+
+void GraphHopcroftTest::testReductionsSameGraphs()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::UndirectedGraph ug2 = getGraph2();
+
+	{
+		// Same simple graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+	}
+
+	{
+		// Same graphs with multi-edges and loops
+		graph::isomorphism::hopcroft::Graph g1(ug2);
+		graph::isomorphism::hopcroft::Graph g2(ug2);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+	}
+}
+
+void GraphHopcroftTest::testReductionsIsomorphicGraphs()
+{
+	graph::UndirectedGraph ug1 = getGraph1();
+	graph::UndirectedGraph ug1_i = getGraph1_I();
+	graph::UndirectedGraph ug2 = getGraph2();
+	graph::UndirectedGraph ug2_i = getGraph2_I();
+
+	{
+		// Isomorphic simple graphs
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1_i);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+	}
+
+	{
+		// Isomorphic graphs with multi-edges and loops
+		graph::isomorphism::hopcroft::Graph g1(ug2);
+		graph::isomorphism::hopcroft::Graph g2(ug2_i);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+	}
+}
+
+void GraphHopcroftTest::testReductionsOther()
+{
+	graph::Node n1("n1");
+	graph::Node n2("n2");
+	graph::Node n3("n3");
+	graph::Node n4("n4");
+	graph::Node n5("n5");
+	graph::Node n6("n6");
+
+	{
+		graph::UndirectedGraph ug1;
+		ug1.addEdge(graph::UndirectedEdge(n1, n5));
+		ug1.addEdge(graph::UndirectedEdge(n1, n6));
+		ug1.addEdge(graph::UndirectedEdge(n1, n2));
+		ug1.addEdge(graph::UndirectedEdge(n1, n3));
+		ug1.addEdge(graph::UndirectedEdge(n2, n3));
+		ug1.addEdge(graph::UndirectedEdge(n2, n4));
+		ug1.addEdge(graph::UndirectedEdge(n2, n5));
+		ug1.addEdge(graph::UndirectedEdge(n2, n6));
+		ug1.addEdge(graph::UndirectedEdge(n3, n4));
+		ug1.addEdge(graph::UndirectedEdge(n3, n5));
+		ug1.addEdge(graph::UndirectedEdge(n4, n5));
+		ug1.addEdge(graph::UndirectedEdge(n5, n6));
+		ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+	}
+
+	{
+		// One node graph
+		graph::UndirectedGraph ug1;
+		ug1.addNode(n1);
+
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug1);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+
+		while (impl.doReductions()) {
+			CPPUNIT_ASSERT_EQUAL(false, impl.m_fail);
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		CPPUNIT_ASSERT_EQUAL(true, impl.testIsomorphism());
+
+		findVertex(&g1, "n1")->vlabel = 1;
+		findVertex(&g2, "n1")->vlabel = 2;
+		CPPUNIT_ASSERT_EQUAL(false, impl.testIsomorphism());
+	}
+
+	{
+		// Random
+		graph::UndirectedGraph ug1 = graph::generate::RandomGraphFactory::generateUndirectedTriConnectedPlanarGraph(30);
+		ug1.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug1));
+		graph::UndirectedGraph ug2 = graph::generate::RandomGraphFactory::generateUndirectedIsomorphicGraph(ug1);
+		ug2.setEmbedding(graph::embedding::HopcroftTarjan::hopcrofttarjan(ug2));
+
+		graph::isomorphism::hopcroft::Graph g1(ug1);
+		graph::isomorphism::hopcroft::Graph g2(ug2);
+		graph::isomorphism::hopcroft::HopcroftImpl impl(&g1, &g2);
+
+		while (impl.doReductions()) {
+			if (impl.m_fail)
+				break;
+			CPPUNIT_ASSERT(checkGraphConsistency(&g1));
+			CPPUNIT_ASSERT(checkGraphConsistency(&g2));
+		}
+		if (impl.m_fail || !impl.testIsomorphism())
+			std::cout << "Random graph fail!" << std::endl;
+	}
+}
diff --git a/alib2algo/test-src/graph/isomorphism/HopcroftTest.h b/alib2algo/test-src/graph/isomorphism/HopcroftTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..ffdf9a483ae044779ee85af4d8e633ae2fa5043c
--- /dev/null
+++ b/alib2algo/test-src/graph/isomorphism/HopcroftTest.h
@@ -0,0 +1,38 @@
+#ifndef HOPCROFT_TEST_H_
+#define HOPCROFT_TEST_H_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class GraphHopcroftTest : public CppUnit::TestFixture
+{
+	CPPUNIT_TEST_SUITE(GraphHopcroftTest);
+	CPPUNIT_TEST(testCircleIsomorphism);
+	CPPUNIT_TEST(testEmbedding);
+	CPPUNIT_TEST(testEdges);
+	CPPUNIT_TEST(testRemoveEdge);
+	CPPUNIT_TEST(testAddInnerEdge);
+	CPPUNIT_TEST(testRemoveDistinguishedEdge);
+	CPPUNIT_TEST(testRemoveFourDegreeExceptionEdge);
+	CPPUNIT_TEST(testAddAndReconnectVertexIntoTriangularFace);
+	CPPUNIT_TEST(testIsomorphismByDefinition);
+	CPPUNIT_TEST(testReductionsSameGraphs);
+	CPPUNIT_TEST(testReductionsIsomorphicGraphs);
+	CPPUNIT_TEST(testReductionsOther);
+	CPPUNIT_TEST_SUITE_END();
+
+public:
+	void testCircleIsomorphism();
+	void testEmbedding();
+	void testEdges();
+	void testRemoveEdge();
+	void testAddInnerEdge();
+	void testRemoveDistinguishedEdge();
+	void testRemoveFourDegreeExceptionEdge();
+	void testAddAndReconnectVertexIntoTriangularFace();
+	void testIsomorphismByDefinition();
+	void testReductionsSameGraphs();
+	void testReductionsIsomorphicGraphs();
+	void testReductionsOther();
+};
+
+#endif // HOPCROFT_TEST_H_
diff --git a/alib2algo/test-src/graph/isomorphism/IsomorphismTest.cpp b/alib2algo/test-src/graph/isomorphism/IsomorphismTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5d5793c2835a69b7b1e7c77980c78d7b0744e8fd
--- /dev/null
+++ b/alib2algo/test-src/graph/isomorphism/IsomorphismTest.cpp
@@ -0,0 +1,239 @@
+#include "IsomorphismTest.h"
+
+#include "graph/isomorphism/Isomorphism.h"
+
+CPPUNIT_TEST_SUITE_NAMED_REGISTRATION(GraphIsomorphismTest, "graph");
+CPPUNIT_TEST_SUITE_REGISTRATION(GraphIsomorphismTest);
+
+void GraphIsomorphismTest::testSameGraphs()
+{
+	graph::Node n1("n1");
+	graph::Node n2("n2");
+	graph::Node n3("n3");
+	graph::Node n4("n4");
+	graph::Node n5("n5");
+
+	// Undirected
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+
+	graph::UndirectedGraph ug2;
+	ug2.addEdge(graph::UndirectedEdge(n1, n2));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge1"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge2"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge3"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n3));
+	ug2.addEdge(graph::UndirectedEdge(n1, n4));
+	ug2.addEdge(graph::UndirectedEdge(n1, n4, "multi-edge4"));
+	ug2.addEdge(graph::UndirectedEdge(n2, n3));
+	ug2.addEdge(graph::UndirectedEdge(n2, n3, "multi-edge5"));
+	ug2.addEdge(graph::UndirectedEdge(n3, n4));
+	ug2.addEdge(graph::UndirectedEdge(n5, n4));
+	ug2.addEdge(graph::UndirectedEdge(n5, n4, "multi-edge6"));
+	ug2.addEdge(graph::UndirectedEdge(n5, n3));
+	ug2.addEdge(graph::UndirectedEdge(n5, n2));
+	ug2.addEdge(graph::UndirectedEdge(n5, n5, "loop1"));
+	ug2.addEdge(graph::UndirectedEdge(n5, n5, "loop2"));
+
+	// Same simple graphs
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(ug1, ug1));
+
+	// Same graphs with multi-edges and loops
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(ug2, ug2));
+
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(ug1, ug2));
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(ug2, ug1));
+
+	// Directed
+	graph::DirectedGraph dg1;
+	dg1.addEdge(graph::DirectedEdge(n1, n2));
+	dg1.addEdge(graph::DirectedEdge(n1, n3));
+	dg1.addEdge(graph::DirectedEdge(n1, n4));
+	dg1.addEdge(graph::DirectedEdge(n2, n3));
+	dg1.addEdge(graph::DirectedEdge(n3, n4));
+	dg1.addEdge(graph::DirectedEdge(n5, n4));
+	dg1.addEdge(graph::DirectedEdge(n5, n3));
+	dg1.addEdge(graph::DirectedEdge(n5, n2));
+
+	graph::DirectedGraph dg2;
+	dg2.addEdge(graph::DirectedEdge(n1, n2));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge1"));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge2"));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge3"));
+	dg2.addEdge(graph::DirectedEdge(n1, n3));
+	dg2.addEdge(graph::DirectedEdge(n1, n4));
+	dg2.addEdge(graph::DirectedEdge(n1, n4, "multi-edge4"));
+	dg2.addEdge(graph::DirectedEdge(n2, n3));
+	dg2.addEdge(graph::DirectedEdge(n2, n3, "multi-edge5"));
+	dg2.addEdge(graph::DirectedEdge(n3, n4));
+	dg2.addEdge(graph::DirectedEdge(n5, n4));
+	dg2.addEdge(graph::DirectedEdge(n5, n4, "multi-edge6"));
+	dg2.addEdge(graph::DirectedEdge(n5, n3));
+	dg2.addEdge(graph::DirectedEdge(n5, n2));
+	dg2.addEdge(graph::DirectedEdge(n5, n5, "loop1"));
+	dg2.addEdge(graph::DirectedEdge(n5, n5, "loop2"));
+
+	// Same simple graphs
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(dg1, dg1));
+
+	// Same graphs with multi-edges and loops
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(dg2, dg2));
+
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(dg1, dg2));
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(dg2, dg1));
+}
+
+void GraphIsomorphismTest::testIsomorphicGraphs()
+{
+	graph::Node n1("n1");
+	graph::Node n2("n2");
+	graph::Node n3("n3");
+	graph::Node n4("n4");
+	graph::Node n5("n5");
+
+	graph::Node n1_("n1_");
+	graph::Node n2_("n2_");
+	graph::Node n3_("n3_");
+	graph::Node n4_("n4_");
+	graph::Node n5_("n5_");
+
+	// Undirected
+	graph::UndirectedGraph ug1;
+	ug1.addEdge(graph::UndirectedEdge(n1, n2));
+	ug1.addEdge(graph::UndirectedEdge(n1, n3));
+	ug1.addEdge(graph::UndirectedEdge(n1, n4));
+	ug1.addEdge(graph::UndirectedEdge(n2, n3));
+	ug1.addEdge(graph::UndirectedEdge(n3, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n4));
+	ug1.addEdge(graph::UndirectedEdge(n5, n3));
+	ug1.addEdge(graph::UndirectedEdge(n5, n2));
+
+	graph::UndirectedGraph ug1_;
+	ug1_.addEdge(graph::UndirectedEdge(n2_, n1_));
+	ug1_.addEdge(graph::UndirectedEdge(n3_, n2_));
+	ug1_.addEdge(graph::UndirectedEdge(n4_, n1_));
+	ug1_.addEdge(graph::UndirectedEdge(n5_, n4_));
+	ug1_.addEdge(graph::UndirectedEdge(n3_, n5_));
+	ug1_.addEdge(graph::UndirectedEdge(n3_, n1_));
+	ug1_.addEdge(graph::UndirectedEdge(n5_, n2_));
+	ug1_.addEdge(graph::UndirectedEdge(n4_, n3_));
+
+	graph::UndirectedGraph ug2;
+	ug2.addEdge(graph::UndirectedEdge(n1, n2));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge1"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge2"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n2, "multi-edge3"));
+	ug2.addEdge(graph::UndirectedEdge(n1, n3));
+	ug2.addEdge(graph::UndirectedEdge(n1, n4));
+	ug2.addEdge(graph::UndirectedEdge(n1, n4, "multi-edge4"));
+	ug2.addEdge(graph::UndirectedEdge(n2, n3));
+	ug2.addEdge(graph::UndirectedEdge(n2, n3, "multi-edge5"));
+	ug2.addEdge(graph::UndirectedEdge(n3, n4));
+	ug2.addEdge(graph::UndirectedEdge(n5, n4));
+	ug2.addEdge(graph::UndirectedEdge(n5, n4, "multi-edge6"));
+	ug2.addEdge(graph::UndirectedEdge(n5, n3));
+	ug2.addEdge(graph::UndirectedEdge(n5, n2));
+	ug2.addEdge(graph::UndirectedEdge(n5, n5, "loop1"));
+	ug2.addEdge(graph::UndirectedEdge(n5, n5, "loop2"));
+
+	graph::UndirectedGraph ug2_;
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n2_));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n2_, "multi-edge1_"));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n2_, "multi-edge2_"));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n2_, "multi-edge3_"));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n3_));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n4_));
+	ug2_.addEdge(graph::UndirectedEdge(n1_, n4_, "multi-edge4_"));
+	ug2_.addEdge(graph::UndirectedEdge(n2_, n3_));
+	ug2_.addEdge(graph::UndirectedEdge(n2_, n3_, "multi-edge5_"));
+	ug2_.addEdge(graph::UndirectedEdge(n3_, n4_));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n4_));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n4_, "multi-edge6_"));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n3_));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n2_));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n5_, "loop1_"));
+	ug2_.addEdge(graph::UndirectedEdge(n5_, n5_, "loop2_"));
+
+	// Isomorphic simple graphs
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(ug1, ug1_));
+
+	// Isomorphic graphs with multi-edges and loops
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(ug2, ug2_));
+
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(ug1, ug2));
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(ug1, ug2_));
+
+	// Directed
+	graph::DirectedGraph dg1;
+	dg1.addEdge(graph::DirectedEdge(n1, n2));
+	dg1.addEdge(graph::DirectedEdge(n1, n3));
+	dg1.addEdge(graph::DirectedEdge(n1, n4));
+	dg1.addEdge(graph::DirectedEdge(n2, n3));
+	dg1.addEdge(graph::DirectedEdge(n3, n4));
+	dg1.addEdge(graph::DirectedEdge(n5, n4));
+	dg1.addEdge(graph::DirectedEdge(n5, n3));
+	dg1.addEdge(graph::DirectedEdge(n5, n2));
+
+	graph::DirectedGraph dg1_;
+	dg1_.addEdge(graph::DirectedEdge(n5_, n4_));
+	dg1_.addEdge(graph::DirectedEdge(n1_, n4_));
+	dg1_.addEdge(graph::DirectedEdge(n1_, n3_));
+	dg1_.addEdge(graph::DirectedEdge(n5_, n2_));
+	dg1_.addEdge(graph::DirectedEdge(n1_, n2_));
+	dg1_.addEdge(graph::DirectedEdge(n3_, n4_));
+	dg1_.addEdge(graph::DirectedEdge(n2_, n3_));
+	dg1_.addEdge(graph::DirectedEdge(n5_, n3_));
+
+	graph::DirectedGraph dg2;
+	dg2.addEdge(graph::DirectedEdge(n1, n2));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge1"));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge2"));
+	dg2.addEdge(graph::DirectedEdge(n1, n2, "multi-edge3"));
+	dg2.addEdge(graph::DirectedEdge(n1, n3));
+	dg2.addEdge(graph::DirectedEdge(n1, n4));
+	dg2.addEdge(graph::DirectedEdge(n1, n4, "multi-edge4"));
+	dg2.addEdge(graph::DirectedEdge(n2, n3));
+	dg2.addEdge(graph::DirectedEdge(n2, n3, "multi-edge5"));
+	dg2.addEdge(graph::DirectedEdge(n3, n4));
+	dg2.addEdge(graph::DirectedEdge(n5, n4));
+	dg2.addEdge(graph::DirectedEdge(n5, n4, "multi-edge6"));
+	dg2.addEdge(graph::DirectedEdge(n5, n3));
+	dg2.addEdge(graph::DirectedEdge(n5, n2));
+	dg2.addEdge(graph::DirectedEdge(n5, n5, "loop1"));
+	dg2.addEdge(graph::DirectedEdge(n5, n5, "loop2"));
+
+	graph::DirectedGraph dg2_;
+	dg2_.addEdge(graph::DirectedEdge(n1_, n2_));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n2_, "multi-edge3_"));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n2_, "multi-edge1_"));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n3_));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n4_, "multi-edge4_"));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n4_));
+	dg2_.addEdge(graph::DirectedEdge(n1_, n2_, "multi-edge2_"));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n4_));
+	dg2_.addEdge(graph::DirectedEdge(n2_, n3_, "multi-edge5_"));
+	dg2_.addEdge(graph::DirectedEdge(n3_, n4_));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n3_));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n5_, "loop1_"));
+	dg2_.addEdge(graph::DirectedEdge(n2_, n3_));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n2_));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n4_, "multi-edge6_"));
+	dg2_.addEdge(graph::DirectedEdge(n5_, n5_, "loop2_"));
+
+	// Isomorphic simple graphs
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(dg1, dg1_));
+
+	// Isomorphic graphs with multi-edges and loops
+	CPPUNIT_ASSERT_EQUAL(true, graph::isomorphism::Isomorphism::isomorphism(dg2, dg2_));
+
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(dg1, dg2));
+	CPPUNIT_ASSERT_EQUAL(false, graph::isomorphism::Isomorphism::isomorphism(dg1, dg2_));
+}
+
diff --git a/alib2algo/test-src/graph/isomorphism/IsomorphismTest.h b/alib2algo/test-src/graph/isomorphism/IsomorphismTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..ff967f2ad5572951765222b72a6f97ff8a07b1a7
--- /dev/null
+++ b/alib2algo/test-src/graph/isomorphism/IsomorphismTest.h
@@ -0,0 +1,18 @@
+#ifndef ISOMORPHISM_TEST_H_
+#define ISOMORPHISM_TEST_H_
+
+#include <cppunit/extensions/HelperMacros.h>
+
+class GraphIsomorphismTest : public CppUnit::TestFixture
+{
+	CPPUNIT_TEST_SUITE(GraphIsomorphismTest);
+	CPPUNIT_TEST(testSameGraphs);
+	CPPUNIT_TEST(testIsomorphicGraphs);
+	CPPUNIT_TEST_SUITE_END();
+
+public:
+	void testSameGraphs();
+	void testIsomorphicGraphs();
+};
+
+#endif // HOPCROFT_TEST_H_
diff --git a/alib2data/src/graph/directed/DirectedEdge.cpp b/alib2data/src/graph/directed/DirectedEdge.cpp
index cfc57db294172bbded2ed09b960c7f7c1174e55b..814355da3f876c00fe4581e86742feaa2a0ecc3a 100644
--- a/alib2data/src/graph/directed/DirectedEdge.cpp
+++ b/alib2data/src/graph/directed/DirectedEdge.cpp
@@ -69,6 +69,10 @@ bool DirectedEdge::operator < (const DirectedEdge& other) const {
 	return this->compare(other) < 0;
 }
 
+bool DirectedEdge::operator > (const DirectedEdge& other) const {
+	return this->compare(other) > 0;
+}
+
 bool DirectedEdge::operator == (const DirectedEdge& other) const {
 	return this->compare(other) == 0;
 }
diff --git a/alib2data/src/graph/directed/DirectedEdge.h b/alib2data/src/graph/directed/DirectedEdge.h
index b8680dbb3d8614c7498c6026a7adc29e13471cbc..8ab4b3f1659efcdb08fbe5d79f111cd07cb15867 100644
--- a/alib2data/src/graph/directed/DirectedEdge.h
+++ b/alib2data/src/graph/directed/DirectedEdge.h
@@ -29,6 +29,7 @@ public:
 	const label::Label &getName() const;
 
 	bool operator < (const DirectedEdge& other) const;
+	bool operator > (const DirectedEdge& other) const;
 	bool operator == (const DirectedEdge& other) const;
 	bool operator != (const DirectedEdge& other) const;
 	int compare(const DirectedEdge &other) const;
diff --git a/alib2data/src/graph/directed/DirectedGraph.cpp b/alib2data/src/graph/directed/DirectedGraph.cpp
index 12ae66f3bab32d633c386a9e9ab71bbbe3a794bf..7822b83c71390c3f48e17f7bbc9fe8add5eeaa91 100644
--- a/alib2data/src/graph/directed/DirectedGraph.cpp
+++ b/alib2data/src/graph/directed/DirectedGraph.cpp
@@ -38,6 +38,7 @@ DirectedGraph::DirectedGraph(const DirectedGraph &other)
 	: representation(other.representation)
 	, nodeValues(other.nodeValues)
 	, edgeValues(other.edgeValues)
+	, embedding(other.embedding)
 {
 	init();
 	impl->copy(other.impl);
@@ -48,6 +49,7 @@ DirectedGraph::DirectedGraph(DirectedGraph &&other) noexcept
 	, impl(other.impl)
 	, nodeValues(std::move(other.nodeValues))
 	, edgeValues(std::move(other.edgeValues))
+	, embedding(std::move(other.embedding))
 {
 	other.impl = 0;
 }
@@ -68,6 +70,7 @@ DirectedGraph &DirectedGraph::operator=(DirectedGraph &&other) noexcept
 	std::swap(this->impl, other.impl);
 	std::swap(this->nodeValues, other.nodeValues);
 	std::swap(this->edgeValues, other.edgeValues);
+	std::swap(this->embedding, other.embedding);
 	return *this;
 }
 
@@ -192,6 +195,16 @@ void DirectedGraph::setEdgeValue(const DirectedEdge &edge, int value)
 	edgeValues.insert({edge, value});
 }
 
+std::unordered_map < Node, std::vector < Node > > DirectedGraph::getEmbedding ( ) const
+{
+	return embedding;
+}
+
+void DirectedGraph::setEmbedding ( const std::unordered_map < Node, std::vector< Node > > embedding )
+{
+	this->embedding = embedding;
+}
+
 void DirectedGraph::operator>>(std::ostream &out) const
 {
 	out << "(DirectedGraph " << static_cast<std::string>(*impl) << ")";
@@ -203,8 +216,8 @@ int DirectedGraph::compare(const DirectedGraph &other) const
 		return static_cast<int>(representation) - static_cast<int>(other.representation);
 	}
 
-	auto first = std::tie(nodeValues, edgeValues);
-	auto second = std::tie(other.nodeValues, other.edgeValues);
+	auto first = std::tie(nodeValues, edgeValues, embedding);
+	auto second = std::tie(other.nodeValues, other.edgeValues, other.embedding);
 
 	std::compare<decltype(first)> comp;
 	int res = comp(first, second);
diff --git a/alib2data/src/graph/directed/DirectedGraph.h b/alib2data/src/graph/directed/DirectedGraph.h
index c828f99ac0546267ea0ea01c4008ae11bcc83a3a..7f1a4efe80851db2d95810e2dc8b4a4e5fc8b16b 100644
--- a/alib2data/src/graph/directed/DirectedGraph.h
+++ b/alib2data/src/graph/directed/DirectedGraph.h
@@ -9,6 +9,7 @@
 #define DIRECTED_GRAPH_H_
 
 #include <unordered_map>
+#include <vector>
 #include <set>
 
 #include "../GraphBase.h"
@@ -62,6 +63,9 @@ public:
 	int getEdgeValue ( const DirectedEdge & edge ) const;
 	void setEdgeValue ( const DirectedEdge & edge, int value );
 
+	std::unordered_map < Node, std::vector < Node > > getEmbedding ( ) const;
+	void setEmbedding ( const std::unordered_map < Node, std::vector< Node > > embedding );
+
 	int compare ( const ObjectBase & other ) const override {
 		if ( std::type_index ( typeid ( * this ) ) == std::type_index ( typeid ( other ) ) ) return this->compare ( ( decltype ( * this ) )other );
 
@@ -88,6 +92,7 @@ private:
 
 	std::unordered_map < Node, int > nodeValues;
 	std::unordered_map < DirectedEdge, int > edgeValues;
+	std::unordered_map < Node, std::vector < Node > > embedding;
 
 	friend class GraphToXMLComposer;
 	friend class GraphToStringComposer;
diff --git a/alib2data/src/graph/undirected/UndirectedEdge.cpp b/alib2data/src/graph/undirected/UndirectedEdge.cpp
index c1d81e73bb2ab4e5ff90d29946d5b13d02bab840..75862596098e61fa0b84d56395b378f2dc152d7a 100644
--- a/alib2data/src/graph/undirected/UndirectedEdge.cpp
+++ b/alib2data/src/graph/undirected/UndirectedEdge.cpp
@@ -74,6 +74,10 @@ bool UndirectedEdge::operator < (const UndirectedEdge& other) const {
 	return this->compare(other) < 0;
 }
 
+bool UndirectedEdge::operator > (const UndirectedEdge& other) const {
+	return this->compare(other) > 0;
+}
+
 bool UndirectedEdge::operator == (const UndirectedEdge& other) const {
 	return this->compare(other) == 0;
 }
diff --git a/alib2data/src/graph/undirected/UndirectedEdge.h b/alib2data/src/graph/undirected/UndirectedEdge.h
index 255f7b30ace1be7e5ebde1b941e4855ca2209e43..5a25b68dbbf3d9c4b5b9e1688e7e4638970c8aa8 100644
--- a/alib2data/src/graph/undirected/UndirectedEdge.h
+++ b/alib2data/src/graph/undirected/UndirectedEdge.h
@@ -29,6 +29,7 @@ public:
 	const label::Label &getName() const;
 
 	bool operator < (const UndirectedEdge& other) const;
+	bool operator > (const UndirectedEdge& other) const;
 	bool operator == (const UndirectedEdge& other) const;
 	bool operator != (const UndirectedEdge& other) const;
 	int compare(const UndirectedEdge &other) const;
diff --git a/alib2data/src/graph/undirected/UndirectedGraph.cpp b/alib2data/src/graph/undirected/UndirectedGraph.cpp
index 3a0ebc79bbcaa88e3e64feb599a4fa8bb9dad5e6..c2f773c37e366fc58f48b53670358e246e2a05d3 100644
--- a/alib2data/src/graph/undirected/UndirectedGraph.cpp
+++ b/alib2data/src/graph/undirected/UndirectedGraph.cpp
@@ -38,6 +38,7 @@ UndirectedGraph::UndirectedGraph(const UndirectedGraph &other)
 	: representation(other.representation)
 	, nodeValues(other.nodeValues)
 	, edgeValues(other.edgeValues)
+	, embedding(other.embedding)
 {
 	init();
 	impl->copy(other.impl);
@@ -48,6 +49,7 @@ UndirectedGraph::UndirectedGraph(UndirectedGraph &&other) noexcept
 	, impl(other.impl)
 	, nodeValues(std::move(other.nodeValues))
 	, edgeValues(std::move(other.edgeValues))
+	, embedding(std::move(other.embedding))
 {
 	other.impl = 0;
 }
@@ -68,6 +70,7 @@ UndirectedGraph &UndirectedGraph::operator=(UndirectedGraph &&other) noexcept
 	std::swap(this->impl, other.impl);
 	std::swap(this->nodeValues, other.nodeValues);
 	std::swap(this->edgeValues, other.edgeValues);
+	std::swap(this->embedding, other.embedding);
 	return *this;
 }
 
@@ -192,6 +195,16 @@ void UndirectedGraph::setEdgeValue(const UndirectedEdge &edge, int value)
 	edgeValues.insert({edge, value});
 }
 
+std::unordered_map < Node, std::vector < Node > > UndirectedGraph::getEmbedding ( ) const
+{
+	return embedding;
+}
+
+void UndirectedGraph::setEmbedding ( const std::unordered_map < Node, std::vector< Node > > embedding )
+{
+	this->embedding = embedding;
+}
+
 void UndirectedGraph::operator>>(std::ostream &out) const
 {
 	out << "(UndirectedGraph " << static_cast<std::string>(*impl) << ")";
@@ -203,8 +216,8 @@ int UndirectedGraph::compare(const UndirectedGraph &other) const
 		return static_cast<int>(representation) - static_cast<int>(other.representation);
 	}
 
-	auto first = std::tie(nodeValues, edgeValues);
-	auto second = std::tie(other.nodeValues, other.edgeValues);
+	auto first = std::tie(nodeValues, edgeValues, embedding);
+	auto second = std::tie(other.nodeValues, other.edgeValues, other.embedding);
 
 	std::compare<decltype(first)> comp;
 	int res = comp(first, second);
diff --git a/alib2data/src/graph/undirected/UndirectedGraph.h b/alib2data/src/graph/undirected/UndirectedGraph.h
index 422b8fdd34c301998393d78bd599b5a2a4b94c6f..5a3d688ba3cf308f10dbb1b7ddb62d0067d7c30b 100644
--- a/alib2data/src/graph/undirected/UndirectedGraph.h
+++ b/alib2data/src/graph/undirected/UndirectedGraph.h
@@ -9,6 +9,7 @@
 #define UNDIRECTED_GRAPH_H_
 
 #include <unordered_map>
+#include <vector>
 #include <set>
 
 #include "../GraphBase.h"
@@ -62,6 +63,9 @@ public:
 	int getEdgeValue ( const UndirectedEdge & edge ) const;
 	void setEdgeValue ( const UndirectedEdge & edge, int value );
 
+	std::unordered_map < Node, std::vector < Node > > getEmbedding ( ) const;
+	void setEmbedding ( const std::unordered_map < Node, std::vector< Node > > embedding );
+
 	int compare ( const ObjectBase & other ) const override {
 		if ( std::type_index ( typeid ( * this ) ) == std::type_index ( typeid ( other ) ) ) return this->compare ( ( decltype ( * this ) )other );
 
@@ -88,6 +92,7 @@ private:
 
 	std::unordered_map < Node, int > nodeValues;
 	std::unordered_map < UndirectedEdge, int > edgeValues;
+	std::unordered_map < Node, std::vector < Node > > embedding;
 
 	friend class GraphToXMLComposer;
 	friend class GraphToStringComposer;