diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1534604..5c07a51 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,7 +3,6 @@ on:
   push:
     branches:
       - main
-      #- dev
     tags: '*'
   pull_request:
 concurrency:
diff --git a/Project.toml b/Project.toml
index 753bd6d..5342acf 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,4 +32,4 @@ SimpleTraits = "0.9"
 SimpleValueGraphs = "0.4"
 SimpleWeightedGraphs = "1.2"
 TensorOperations = "3.2"
-julia = "1.8"
\ No newline at end of file
+julia = "1.8"
diff --git a/README.md b/README.md
index 13cbc4d..f9758b7 100644
--- a/README.md
+++ b/README.md
@@ -39,7 +39,6 @@ In the package documentation you can find a [tutorial](https://juliagraphs.org/M
 
 ## Future Developments 
 
-- [ ] [Implement faster graph realization algorithms](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/32);
 - [ ] [Implement more general configuration models / graph generators](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/33);
 - [ ] [Implement graph of layers](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/34);
 - [ ] [Implement projected monoplex and overlay graphs](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/35);
diff --git a/docs/src/index.md b/docs/src/index.md
index edf3117..372677d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -852,7 +852,6 @@ Read a complete list of analytical methods exclusive to multilayer graphs in the
 
 ### Future Developments
 
-- [Implement faster graph realization algorithms](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/32);
 - [Implement more general configuration models / graph generators](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/33);
 - [Implement graph of layers](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/34);
 - [Implement projected monoplex and overlay graphs](https://github.com/JuliaGraphs/MultilayerGraphs.jl/issues/35);
diff --git a/mains/abstractmultilayerugraph_development.jl b/mains/abstractmultilayerugraph_development.jl
new file mode 100644
index 0000000..66f57eb
--- /dev/null
+++ b/mains/abstractmultilayerugraph_development.jl
@@ -0,0 +1,381 @@
+using Revise, Test, Logging, LoggingExtras, Beep
+using DataStructures
+using StatsBase, Distributions #, Bijections, SparseArrays
+using Graphs, SimpleWeightedGraphs, MetaGraphs, SimpleValueGraphs
+using Agents
+using MultilayerGraphs
+
+
+degree_sequence = degree(SimpleGraph(10,20))
+
+havel_hakimi_graph_generator_unsafe(degree_sequence)
+
+graph = SimpleGraph(length(degree_sequence))
+v_ds = OrderedDict(i => deg for (i,deg) in enumerate(degree_sequence))
+
+
+while(any(values(v_ds) .!= 0))
+    v_ds = OrderedDict(sort(collect(v_ds), by = last , rev = true))
+    S,s = popfirst!(v_ds)
+    @debug "" S, s
+    all(collect(values(v_ds))[1:s] .> 0) || throw(ErrorException("The provided `degree_sequence` is not graphical"))
+    # println("here")
+    for v in collect(keys(v_ds))[1:s]
+        add_edge!(graph, S, v)
+        @debug "", v, length(v_ds), collect(keys(v_ds))
+        v_ds[v] -= 1
+        println("here")
+    end
+end
+
+degree(graph) .== degree_sequence
+
+# cd("dev/MultilayerGraphs")
+# cd("dev/MultilayerGraphs/test")
+
+# Create the logger ushc that it prints debug messages only coming from our package. Taken from https://julialogging.github.io/how-to/filter-messages/#How-to-filter-messages
+logger = EarlyFilteredLogger(ConsoleLogger(stderr, Logging.Debug)) do args
+    r =  args._module === MultilayerGraphs || args._module === Main
+    return r
+end
+global_logger(logger)
+
+
+const vertextype   = Int64
+const _weighttype  = Float64 
+const min_vertices = 5
+const max_vertices = 7
+const min_edges    = 1
+const max_edges    = max_vertices*(max_vertices-1)
+
+const multilayer_nodes = [Node("node_$i") for i in 1:max_vertices]
+const mvs_layers = MV.(multilayer_nodes)
+# Multilayer vertices with metadata (to be used with MetaGraphs, SimpleValueGraphs, etc)
+const mvs_metadata  = [MV(node, ("I'm node $(node.id)",)) for node in multilayer_nodes]
+
+function _get_srcmv_dstmv_layer(layer::Layer)
+
+    mvs = MultilayerGraphs.get_bare_mv.(collect(mv_vertices(layer)))
+
+    src_mv = nothing    
+    _collection = []
+
+
+    while isempty(_collection)
+        src_mv = rand(mvs)
+        _collection = setdiff(Set(mvs), Set(vcat(MultilayerGraphs.get_bare_mv.(mv_outneighbors(layer, src_mv)), src_mv ) ) )  
+    end
+
+    dst_mv = MultilayerGraphs.get_bare_mv(rand(_collection))
+
+    return mvs, src_mv, dst_mv
+end
+
+
+_nv = rand(min_vertices:max_vertices)
+_ne =  rand(_nv:(_nv*(_nv-1)) ÷ 2 )
+layer_sg = Layer(:layer_sg, sample(mvs_layers, _nv, replace = false), _ne, SimpleGraph{vertextype}(), _weighttype)
+
+
+_nv = rand(min_vertices:max_vertices)
+_ne =  rand(_nv:(_nv*(_nv-1)) ÷ 2 )
+layer_swg = Layer(:layer_swg, sample(mvs_layers, _nv, replace = false), _ne, SimpleWeightedGraph{vertextype, _weighttype}(), _weighttype;  default_edge_weight = (src,dst) -> rand())
+
+
+_nv = rand(min_vertices:max_vertices)
+_ne =  rand(_nv:(_nv*(_nv-1)) ÷ 2 )
+
+layer_mg = Layer(:layer_mg, sample(mvs_metadata, _nv, replace = false), _ne, MetaGraph{vertextype, _weighttype}(), _weighttype; default_edge_metadata = (src,dst) -> (from_to = "from_$(src)_to_$(dst)",))
+
+
+_nv = rand(min_vertices:max_vertices)
+_ne =  rand(_nv:(_nv*(_nv-1)) ÷ 2 )
+layer_vg = Layer(:layer_vg, 
+                 sample(mvs_metadata, _nv, replace = false), 
+                 _ne,
+                 MultilayerGraphs.ValGraph(SimpleGraph{vertextype}();edgeval_types=(Float64, String, ),
+                                            edgeval_init=(s, d) -> (s+d, "hi"),
+                                            vertexval_types=(String,),
+                                            vertexval_init=v -> ("$v",),),
+                _weighttype;
+                default_edge_metadata = (src,dst) -> (rand(), "from_$(src)_to_$(dst)"),
+                default_vertex_metadata = mv -> ("This metadata have been generated via the default_vertex_metadata method",))
+
+
+
+all_layers_u = [layer_sg, layer_swg, layer_mg, layer_vg];
+
+
+_nv = nv(layer_sg) + nv(layer_swg)
+_ne = rand(_nv:(_nv*(_nv-1))÷ 2)
+interlayer_sg_swg = Interlayer(layer_sg, layer_swg, _ne,  SimpleGraph{vertextype}())
+
+_nv =  nv(layer_swg) + nv(layer_mg)
+_ne = rand(_nv:(_nv*(_nv-1))÷ 2)
+interlayer_swg_mg = Interlayer(layer_swg, layer_mg, _ne, SimpleWeightedGraph{vertextype, _weighttype}(); default_edge_weight = (x,y) -> rand() ) 
+
+_nv =  nv(layer_mg) + nv(layer_vg)
+_ne = rand(_nv:(_nv*(_nv-1))÷ 2)
+interlayer_mg_vg = Interlayer(layer_mg, layer_vg, _ne, MetaGraph{vertextype, _weighttype}(); default_edge_metadata = (x,y) -> (mymetadata = rand(),), transfer_vertex_metadata = true )  #SimpleWeightedGraph{Int64, Float64}()
+
+interlayer_multiplex_sg_mg = multiplex_interlayer(layer_sg, layer_mg, ValGraph(SimpleGraph{vertextype}(); edgeval_types=(from_to = String,), edgeval_init=(s, d) -> (from_to = "from_$(s)_to_$(d)")); default_edge_metadata = (x,y) -> (from_to = "from_$(src)_to_$(dst)",)) 
+
+interlayer_empty_sg_vg = empty_interlayer(layer_sg, layer_vg, SimpleGraph{vertextype}())
+
+
+
+all_interlayers_u = [interlayer_sg_swg, interlayer_swg_mg, interlayer_mg_vg, interlayer_multiplex_sg_mg, interlayer_empty_sg_vg]
+
+
+multilayergraph = MultilayerGraph(all_layers_u, all_interlayers_u) #, all_interlayers_u
+
+
+# Test random multilayergraph (configuration model)
+
+# We need to wait for this PR: https://github.com/JuliaGraphs/SimpleWeightedGraphs.jl/pull/14 top be merged before being able to use SimpleWeightedGraphs in the configuration model. Also we cannot use AbstractValgraphs until we implement default weight and metadata functions for layers and interlayers to be kept into the descriptor.
+
+layers_to_be_emptied =  deepcopy([layer for layer in all_layers_u if !(layer.graph isa SimpleWeightedGraphs.AbstractSimpleWeightedGraph)])
+
+layers_names_to_be_emptied = name.(layers_to_be_emptied)
+
+interlayers_to_be_emptied =  deepcopy([interlayer for interlayer in all_interlayers_u if all(in.(interlayer.layers_names, Ref(layers_names_to_be_emptied))) && !(interlayer.graph isa SimpleWeightedGraphs.AbstractSimpleWeightedGraph) ])
+
+
+for layer in layers_to_be_emptied
+    for edge in edges(layer)
+        rem_edge!(layer, edge)
+    end
+end
+
+
+for interlayer in interlayers_to_be_emptied
+    for edge in edges(interlayer)
+        rem_edge!(interlayer, edge)
+    end
+end
+
+
+@test all(ne.(layers_to_be_emptied) .== 0)
+
+@test all(ne.(interlayers_to_be_emptied) .== 0)
+
+
+# Instantiate configuration-model multilayergraph
+
+# configuration_multilayergraph = MultilayerGraph(layers_to_be_emptied, interlayers_to_be_emptied, truncated(Normal(10), 0.0, 20.0) );
+
+# histogram(degree(configuration_multilayergraph), bins = 20)
+
+
+# Test layers and interlayers
+@test nl(multilayergraph) == length(all_layers_u)
+_nl = nl(multilayergraph) 
+@test nIn(multilayergraph) == _nl*(_nl-1) / 2
+
+# Test get_interlayer
+for (layer_1, layer_2) in Iterators.product(multilayergraph.layers_names, multilayergraph.layers_names)
+    if layer_1 != layer_2
+        interlayer = get_interlayer(multilayergraph,  layer_1, layer_2)
+        @test nv(interlayer) >= 0
+    end
+end
+
+# Test add_layer! and rem_layer!
+const n_vertices_missing = rand(min_vertices:max_vertices)
+const n_edges_missing_u =  rand(n_vertices_missing:(n_vertices_missing*(n_vertices_missing-1))÷ 2)
+missing_layer_u = Layer(:missing_layer_u, 
+sample(mvs_metadata, n_vertices_missing, replace = false), 
+n_edges_missing_u,
+MultilayerGraphs.ValGraph(SimpleGraph{vertextype}();edgeval_types=(Float64, String, ),
+                           edgeval_init=(s, d) -> (s+d, "missing vertex $(s+d)"),
+                           vertexval_types=(String,),
+                           vertexval_init=v -> ("$(v^2)",),),
+_weighttype;
+default_edge_metadata = (src,dst) -> (rand(), "missing edge from_$(src)_to_$(dst)",)
+) # SimpleGraph{vertextype}()
+
+@test !has_layer(multilayergraph, missing_layer_u.name)
+@test add_layer!(multilayergraph, missing_layer_u)
+@test has_layer(multilayergraph, missing_layer_u.name)
+@test rem_layer!(multilayergraph, missing_layer_u.name)
+@test !has_layer(multilayergraph, missing_layer_u.name)
+
+# Test nodes
+@inferred(nodes(multilayergraph))
+
+@inferred(nn(multilayergraph))
+
+for node in nodes(multilayergraph)
+    @test has_node(multilayergraph, node)
+end
+
+## Test MultilayerGraphs.add_node! and MultilayerGraphs.rem_node!
+new_node = Node("new_node")
+nv_prev = nv(multilayergraph)
+ne_prev = ne(multilayergraph)
+@test !has_node(multilayergraph, new_node)
+@test MultilayerGraphs.add_node!(multilayergraph, new_node)
+@test has_node(multilayergraph, new_node)
+@test MultilayerGraphs.rem_node!(multilayergraph, new_node)
+@test !has_node(multilayergraph, new_node)
+### Test that nothing changed
+@test nv_prev == nv(multilayergraph)
+@test ne_prev == ne(multilayergraph)
+
+# Test vertices
+@test eltype(multilayergraph) == Int64
+nv(multilayergraph)
+@test  length(multilayergraph.fadjlist) == length(vertices(multilayergraph)) # nv_withmissing(multilayergraph) ==
+
+
+## Test that all multilayer vertices are present
+for mv in vcat(mv_vertices.(all_layers_u)...)
+    @test has_vertex(multilayergraph, mv)
+end
+
+for mv in mv_vertices(multilayergraph)
+    # if !(mv isa MissingVertex)
+        mv_inneighbors(multilayergraph, mv)
+        mv_outneighbors(multilayergraph, mv)
+    # end
+end
+
+## Test set_edge_weight!
+_, rand_mv_1_weight, rand_mv_2_weight = _get_srcmv_dstmv_layer(layer_swg)
+_weight = 3.14
+@test !has_edge(multilayergraph, rand_mv_1_weight, rand_mv_2_weight)
+@test add_edge!(multilayergraph, rand_mv_1_weight, rand_mv_2_weight, weight = _weight)
+@test has_edge(multilayergraph, rand_mv_1_weight, rand_mv_2_weight)
+wgt = weight_tensor(multilayergraph)
+@test wgt[rand_mv_1_weight, rand_mv_2_weight] == wgt[rand_mv_2_weight, rand_mv_1_weight] == get_weight(multilayergraph, rand_mv_1_weight, rand_mv_2_weight) == get_weight(multilayergraph, rand_mv_2_weight, rand_mv_1_weight)  == _weight
+@test set_weight!(multilayergraph , rand_mv_1_weight, rand_mv_2_weight, _weight + 1)
+wgt = weight_tensor(multilayergraph)
+@test wgt[rand_mv_1_weight, rand_mv_2_weight] == wgt[rand_mv_2_weight, rand_mv_1_weight] == get_weight(multilayergraph, rand_mv_1_weight, rand_mv_2_weight) == get_weight(multilayergraph, rand_mv_2_weight, rand_mv_1_weight)== _weight + 1
+
+## Test set_metadata!
+_, rand_mv_1_meta, rand_mv_2_meta = _get_srcmv_dstmv_layer(layer_mg)
+### On vertices
+@test set_metadata!(multilayergraph, rand_mv_1_meta, (meta = "new_metadata",))
+@test get_metadata(multilayergraph, rand_mv_1_meta).meta == "new_metadata"
+### On edges
+@test !has_edge(multilayergraph, rand_mv_1_meta, rand_mv_2_meta)
+@test add_edge!(multilayergraph, rand_mv_1_meta, rand_mv_2_meta, metadata =  (meta = "hello",))
+@test has_edge(multilayergraph, rand_mv_1_meta, rand_mv_2_meta)
+mt = metadata_tensor(multilayergraph)
+@test mt[rand_mv_1_meta, rand_mv_2_meta].meta == mt[rand_mv_2_meta, rand_mv_1_meta].meta == get_metadata(multilayergraph, rand_mv_1_meta, rand_mv_2_meta).meta == get_metadata(multilayergraph, rand_mv_2_meta, rand_mv_1_meta).meta == "hello"
+_metadata = (meta = "byebye",)
+@test set_metadata!(multilayergraph , rand_mv_1_meta, rand_mv_2_meta, _metadata)
+mt = metadata_tensor(multilayergraph)
+@test mt[rand_mv_1_meta, rand_mv_2_meta].meta == mt[rand_mv_2_meta, rand_mv_1_meta].meta == get_metadata(multilayergraph, rand_mv_1_meta, rand_mv_2_meta).meta == get_metadata(multilayergraph, rand_mv_2_meta, rand_mv_1_meta).meta == "byebye"
+
+
+
+
+
+
+
+
+# Test edges
+ne(multilayergraph)
+
+## Test that all edges are present
+for edge in vcat(collect.(edges.(all_layers_u))..., collect.(edges.(all_interlayers_u))... )
+    @test has_edge(multilayergraph, edge)
+end
+
+@test MultilayerGraphs.weighttype(multilayergraph) == Float64
+@test edgetype(multilayergraph) == MultilayerEdge{Float64}
+
+# Test set_edge_weight!
+
+
+## Test add_edge! and rem_edge!
+
+
+# Test Graphs.jl's extra overrides
+@test all(indegree(multilayergraph)  .== degree(multilayergraph)) #.+ outdegree(multilayergraph)
+
+@inferred(mean_degree(multilayergraph))
+@inferred(degree_second_moment(multilayergraph))
+@inferred(degree_variance(multilayergraph))
+
+
+
+
+# Test multilayer-specific methods
+@test all(get_supra_weight_matrix_from_weight_tensor(weight_tensor(multilayergraph).array) .== supra_weight_matrix(multilayergraph).array)
+@test all(get_weight_tensor_from_supra_weight_matrix(multilayergraph, supra_weight_matrix(multilayergraph).array) .==  weight_tensor(multilayergraph).array)
+@test_broken multilayer_global_clustering_coefficient(multilayergraph) .== global_clustering_coefficient(multilayergraph)
+
+overlaygraph = get_overlay_monoplex_graph(multilayergraph)
+@test_broken global_clustering_coefficient(overlaygraph) .== overlay_clustering_coefficient(multilayergraph)
+
+@test multilayer_weighted_global_clustering_coefficient(multilayergraph, [1 / 3, 1 / 3, 1 / 3]) .≈ multilayer_global_clustering_coefficient(multilayergraph)
+
+eig_centr_u, errs_u = eigenvector_centrality(multilayergraph; norm="n", tol=1e-3)
+
+modularity(multilayergraph, rand([1, 2, 3, 4], length(nodes(multilayergraph)), length(multilayergraph.layers)), )
+
+get_graph_of_layers(multilayergraph)
+
+
+wgt = weight_tensor(multilayergraph)
+sam = supra_weight_matrix(multilayergraph)
+for edge in collect(edges(multilayergraph.layer_swg))
+    @debug "" src(edge) dst(edge) edge
+    @test wgt[src(edge), dst(edge)] == wgt[ dst(edge), src(edge)] == MultilayerGraphs.weight(edge)
+    @test sam[src(edge), dst(edge)] ==  sam[ dst(edge), src(edge)] ==  MultilayerGraphs.weight(edge)
+end
+
+
+
+
+
+
+# Test that, given a 1-dimensional multilayergraph, we obtain the same metrics as we would by using Graphs.jl's utilities on the one and only layer
+## unweighted and weighted case
+for layer in all_layers_u
+
+    if !(layer.graph isa SimpleValueGraphs.AbstractValGraph)
+        monolayergraph = MultilayerGraph([layer])
+
+        @test length(edges(monolayergraph)) == length(edges(layer.graph))
+
+        @test eltype(monolayergraph) == eltype(layer.graph)
+
+        @test ne(monolayergraph) == ne(layer.graph)
+
+        @test length(nodes(monolayergraph)) == nv(layer.graph)
+
+        @test nv(monolayergraph) .== nv(layer.graph)
+
+        @test all(inneighbors.(Ref(monolayergraph), vertices(monolayergraph)) .== inneighbors.(Ref(layer.graph), vertices(layer.graph)))
+
+
+        @test all(indegree(monolayergraph) .== indegree(layer.graph))
+
+        @test all(outdegree(monolayergraph) .== outdegree(layer.graph))
+
+        @test all(degree(monolayergraph) .== degree(layer.graph))
+
+        @test_broken vec(eigenvector_centrality(monolayergraph; norm="n", tol=1e-3)[1]) ==
+            eigenvector_centrality(layer.graph)
+
+
+        tests = Bool[]
+        for i in 1:5
+            clustering = rand(
+                [1, 2, 3], length(nodes(monolayergraph)), length(monolayergraph.layers)
+            )
+            push!(
+                tests,
+                modularity(monolayergraph, clustering) == modularity(layer.graph, vec(clustering)),
+            )
+        end
+        @test_broken all(tests)
+
+        for edge in edges(layer)
+            @test has_edge(monolayergraph, edge)
+        end
+    end
+end
\ No newline at end of file
diff --git a/src/MultilayerGraphs.jl b/src/MultilayerGraphs.jl
index 66f4004..4e6a589 100644
--- a/src/MultilayerGraphs.jl
+++ b/src/MultilayerGraphs.jl
@@ -146,7 +146,9 @@ export
     δ_1,
     δ_2,
     δ_3,
-    δ_Ω
+    δ_Ω,
+    havel_hakimi_graph_generator,
+    kleitman_wang_graph_generator
     # tensorfacoriazations.jl
 
 using Base, InteractiveUtils, IterTools, SimpleTraits, Bijections
diff --git a/src/multilayerdigraph.jl b/src/multilayerdigraph.jl
index b1513ec..fd61711 100644
--- a/src/multilayerdigraph.jl
+++ b/src/multilayerdigraph.jl
@@ -210,7 +210,10 @@ function MultilayerDiGraph(
         
     end
     
-    edge_list = _random_directed_configuration(_multilayerdigraph, indegree_sequence, outdegree_sequence, allow_self_loops)
+    # edge_list = _random_directed_configuration(_multilayerdigraph, indegree_sequence, outdegree_sequence, allow_self_loops)
+    equivalent_graph = kleitman_wang_graph_generator(indegree_sequence, outdegree_sequence)
+
+    edge_list = [ME(_multilayerdigraph.v_V_associations[src(edge)], _multilayerdigraph.v_V_associations[dst(edge)]) for edge in edges(equivalent_graph) ]
 
     for edge in edge_list
         add_edge!(_multilayerdigraph, edge)
diff --git a/src/multilayergraph.jl b/src/multilayergraph.jl
index e1254ea..bc492b8 100644
--- a/src/multilayergraph.jl
+++ b/src/multilayergraph.jl
@@ -202,7 +202,10 @@ function MultilayerGraph(empty_multilayergraph::MultilayerGraph{T,U}, degree_seq
         isgraphical(degree_sequence) || throw(ArgumentError("degree_sequence must be graphical.")) 
     end
 
-    edge_list = _random_undirected_configuration(_multilayergraph, degree_sequence, allow_self_loops)
+    # edge_list = _random_undirected_configuration(_multilayergraph, degree_sequence, allow_self_loops)
+    equivalent_graph = havel_hakimi_graph_generator(degree_sequence)
+
+    edge_list = [ME(empty_multilayergraph.v_V_associations[src(edge)], empty_multilayergraph.v_V_associations[dst(edge)]) for edge in edges(equivalent_graph) ]
 
     for edge in edge_list
         add_edge!(_multilayergraph, edge)
diff --git a/src/utilities.jl b/src/utilities.jl
index 356ba14..2b93b2c 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -1,10 +1,10 @@
 # Base extensions
-"""
+#= """
     getindex(od::O, key::Int64) where {O <: OrderedDict}
 
 Return the `key`th pair of `od`, following `od`'s order.
 """
-Base.getindex(od::O, key::Int64) where {O<:OrderedDict} = collect(values(od))[key]
+Base.getindex(od::O, key::Int64) where {O<:OrderedDict} = collect(values(od))[key] =#
 
 # Non package-specific utilities
 """
@@ -426,7 +426,7 @@ function isdigraphical(
 end
 
 
-"""
+#= """
     _random_undirected_configuration(empty_mg::M, degree_sequence::Vector{ <: Integer}) where {T,U,M <: MultilayerGraph{T,U}}
 
 Internal function. Returns a `MultilayerEdge` list compatible with `empty_mg`, using a relatively inefficient algorithm.
@@ -548,7 +548,7 @@ function _random_directed_configuration(empty_mg::M, indegree_sequence::Vector{
             success = length(mvs_indegree_dict) == 0 && length(mvs_outdegree_dict) == 0
         end
     return edge_list
-end
+end =#
 
 
 """
@@ -558,10 +558,8 @@ Internal function. Converts `cart_index` to an integer index such that it corres
 """
 cartIndexTovecIndex(cart_index::Union{NTuple{N, Integer},CartesianIndex}, tensor_size::NTuple{N, <: Integer} ) where N = cart_index[1] + sum( collect(Tuple(cart_index)[2:end] .- 1) .* cumprod(tensor_size[1:end-1]))
 
-
-# TODO:
 """
-    havel_hakimi_(empty_graph::SimpleGraph, degree_sequence::Vector{<:Integer})
+    havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer})
 
 Returns a simple graph with a given finite degree sequence of non-negative integers generated via the Havel-Hakimi algorithm which works as follows:
 1. successively connect the node of highest degree to other nodes of highest degree;
@@ -569,21 +567,123 @@ Returns a simple graph with a given finite degree sequence of non-negative integ
 3. repeat the procedure.
 
 ## References
-1. [Hakimi (1962)](https://doi.org/10.1137/0110037)
-2. [Kleitman and Wang (1973)](https://doi.org/10.1016/0012-365X(73)90037-X)
-"""
-function havel_hakimi_(empty_graph::SimpleGraph, degree_sequence::Vector{<:Integer}) # Please think about a decent name!
-    # Check whether the given degree sequence contains only non-negative integers
-    !any(degree -> degree < 0, degree_sequence) || throw(ArgumentError("The degree sequence (degree_sequence) is invalid: it must contain non-negative integers only."))
-    # Check whether the given degree sequence is compatible with the given multilayer graph
-    nv(empty_graph) == length(degree_sequence) || throw(ArgumentError("The degree sequence (degree_sequence) and the multilayer graph (empty_mg) are incompatible: the length of the degree sequence doesn't coincide with the number of vertices."))
-    # Check whether the given degree sequence is graphical
-    isgraphical(degree_sequence) || throw(ArgumentError("The degree sequence (degree_sequence) is invalid: it must be graphical (i.e. realizable in a simple graph)."))
-    # Check whether the given multilayer graph is undirected
-    !is_directed(empty_mg) || throw(ArgumentError("The multilayer graph (empty_mg) is invalid: it must be undirected."))
-    # Get all the multilayer vertices from the empty multilayer graph
-    ##mvs = mv_vertices(empty_mg)
-    # Get the length of the degree sequence
-    ##n = length(degree_sequence)
-    ##mvs_degree_dict = Dict(mv => deg for (mv, deg) in zip(mvs, degree_sequence) if deg > 0)
+1. [Hakimi (1962)](https://doi.org/10.1137/0110037);
+2. [Wikipedia](https://en.wikipedia.org/wiki/Havel%E2%80%93Hakimi_algorithm).
+"""
+function havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer})
+    # Check whether the degree sequence has only non-negative values
+    all(degree_sequence .>= 0) || throw(ArgumentError("The degree sequence must contain non-negative integers only."))
+    # Instantiate an empty simple graph 
+    graph = SimpleGraph(length(degree_sequence))
+    # Create a (vertex, degree) ordered dictionary
+    vertices_degrees_dict = OrderedDict(vertex => degree for (vertex, degree) in enumerate(degree_sequence))
+    # Havel-Hakimi algorithm  
+    while(any(values(vertices_degrees_dict) .!= 0))
+        # Sort the new sequence in non-increasing order
+        vertices_degrees_dict = OrderedDict(sort(collect(vertices_degrees_dict), by = last , rev = true))
+        # Remove the first vertex and distribute its stabs
+        max_vertex, max_degree = popfirst!(vertices_degrees_dict)
+        # Check whether the new sequence has only positive values
+        all(collect(values(vertices_degrees_dict))[1:max_degree] .> 0) || throw(ErrorException("The degree sequence is not graphical."))
+        # Connect the node of highest degree to other nodes of highest degree 
+        for vertex in collect(keys(vertices_degrees_dict))[1:max_degree]
+            add_edge!(graph, max_vertex, vertex)
+            vertices_degrees_dict[vertex] -= 1
+        end
+    end
+    # Return the simple graph
+    return graph
+end
+
+"""
+    lexicographical_order_lt(A::Vector{T}, B::Vector{T}) where T
+
+The less than (lt) function that implements lexicographical order.
+
+See [Wikipedia](https://en.wikipedia.org/wiki/Lexicographic_order).
+"""
+function lexicographical_order_lt(A::Union{Vector{T},NTuple{N,T}}, B::Union{Vector{T}, NTuple{M,T}}) where {N,M,T}
+    A_dc = deepcopy(A)
+    B_dc = deepcopy(B)
+    diff_length = length(A) - length(B)
+    if diff_length >= 0
+        A_dc = vcat(A_dc, repeat([-Inf], diff_length))
+    else
+        B_dc = vcat(B_dc, repeat([-Inf], abs(diff_length)))
+    end
+    for (a,b) in zip(A_dc, B_dc)
+        if a != b
+            a < b
+        end
+    end
 end
+
+"""
+lexicographical_order_ntuple(A::NTuple{N,T}, B::NTuple{M,T}) where {N,T}
+
+The less than (lt) function that implements lexicographical order for `NTuple` of equal length.
+
+See [Wikipedia](https://en.wikipedia.org/wiki/Lexicographic_order).
+"""
+function lexicographical_order_ntuple(A::NTuple{N,T}, B::NTuple{N,T}) where {N,T}
+
+    for (a,b) in zip(A, B)
+        if a != b
+            return a < b
+        end
+    end
+
+    return false
+
+end
+
+"""
+    kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_sequence::AbstractVector{<:Integer})
+
+Returns a simple directed graph with given finite in-degree and out-degree sequences of non-negative integers generated via the Kleitman-Wang algorithm, that works like follows:
+1. Sort the indegree-outdegree pairs in lexicographical order;
+2. Select a pair that has strictly positive outdegree, say the i-th pairs that has outdegree = b_i;
+3. Subtract 1 to the first b_i highest indegrees (the i-th being excluded), and set b_i to 0;
+4. Repeat from 1. until all indegree-outdegree pairs are of the form (0.0).
+
+## References
+- [Wikipedia](https://en.wikipedia.org/wiki/Kleitman%E2%80%93Wang_algorithms)
+- [Kleitman and Wang (1973)](https://doi.org/10.1016/0012-365X(73)90037-X)
+"""
+function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer}, outdegree_sequence::AbstractVector{<:Integer})
+
+    length(indegree_sequence) == length(outdegree_sequence) || throw(ArgumentError("The provided `indegree_sequence` and `outdegree_sequence` must be of the dame length."))
+    # Check whether the indegree_sequence and outdegree_sequence have only non-negative values
+    all(indegree_sequence .>= 0) || throw(ArgumentError("The `indegree_sequence` sequence must contain non-negative integers only."))
+    all(outdegree_sequence .>= 0) || throw(ArgumentError("The `outdegree_sequence` sequence must contain non-negative integers only."))
+
+    # Instantiate an empty simple graph 
+    graph = SimpleDiGraph(length(indegree_sequence))
+    # Create a (vertex, degree) ordered dictionary
+    S = zip(deepcopy(indegree_sequence), deepcopy(outdegree_sequence))
+    vertices_degrees_dict = OrderedDict(i => tup for (i,tup) in enumerate(S))
+    # Kleitman-Wang algorithm  
+    while(any(Iterators.flatten(values(vertices_degrees_dict)) .!= 0 ))
+        # Sort the new sequence in non-increasing lexicographical order
+        vertices_degrees_dict = OrderedDict(sort(collect(vertices_degrees_dict), by = last, lt = lexicographical_order_ntuple , rev = true))
+        # Find a vertex with positive outdegree,a nd temporarily remove it from `vertices_degrees_dict`
+        i, (a_i, b_i) = 0,(0,0)
+        for (_i, (_a_i, _b_i)) in collect(deepcopy(vertices_degrees_dict))
+            if _b_i != 0
+                i, a_i, b_i = (_i,_a_i, _b_i)
+                delete!(vertices_degrees_dict, _i)
+                break
+            end
+        end        
+        # Connect the vertex found above to other nodes of highest degree 
+        for (v,degs) in collect(vertices_degrees_dict)[1:b_i]
+            add_edge!(graph, i, v)
+            vertices_degrees_dict[v] = (degs[1]-1, degs[2])
+        end
+        # Check whether the new sequence has only positive values
+        all(collect(Iterators.flatten(collect(values(vertices_degrees_dict))))[1:b_i] .>= 0) || throw(ErrorException("The in-degree and out-degree sequences are not digraphical."))
+        # Reinsert the vertex, with zero outdegree
+        vertices_degrees_dict[i] = (a_i, 0)
+    end
+    return graph
+end
\ No newline at end of file
diff --git a/test/abstractmultilayerdigraph.jl b/test/abstractmultilayerdigraph.jl
index 2ac659e..983a9df 100644
--- a/test/abstractmultilayerdigraph.jl
+++ b/test/abstractmultilayerdigraph.jl
@@ -21,7 +21,7 @@ end
 @test all(ne.(interlayers_to_be_emptied) .== 0)
 
 # Instantiate configuration-model multilayerdigraph
-# configuration_multilayerDigraph = MultilayerDiGraph(layers_to_be_emptied, interlayers_to_be_emptied, truncated(Normal(10), 0.0, 20.0), truncated(Normal(11), 0.0, 22.0));
+configuration_multilayerDigraph = MultilayerDiGraph(layers_to_be_emptied, interlayers_to_be_emptied, truncated(Normal(10), 0.0, 20.0), truncated(Normal(11), 0.0, 22.0));
 
 # Test get_interlayer
 for (layer_1, layer_2) in Iterators.product(multilayerdigraph.layers_names, multilayerdigraph.layers_names)
diff --git a/test/abstractmultilayerugraph.jl b/test/abstractmultilayerugraph.jl
index ca16b5b..972e236 100644
--- a/test/abstractmultilayerugraph.jl
+++ b/test/abstractmultilayerugraph.jl
@@ -21,7 +21,7 @@ end
 @test all(ne.(interlayers_to_be_emptied) .== 0)
 
 # Instantiate configuration-model multilayergraph
-# configuration_multilayergraph = MultilayerGraph(layers_to_be_emptied, interlayers_to_be_emptied, truncated(Normal(10), 0.0, 20.0));
+configuration_multilayergraph = MultilayerGraph(layers_to_be_emptied, interlayers_to_be_emptied, truncated(Normal(10), 0.0, 20.0));
 
 # Test get_interlayer
 for (layer_1, layer_2) in Iterators.product(multilayergraph.layers_names, multilayergraph.layers_names)