From 572f560986569664c31ba091c864ce08edbbcb92 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Wed, 14 Dec 2022 22:12:37 +0100
Subject: [PATCH 01/17] Implement havel-akimi (currently fails)

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/MultilayerGraphs.jl |  3 ++-
 src/utilities.jl        | 44 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 46 insertions(+), 1 deletion(-)

diff --git a/src/MultilayerGraphs.jl b/src/MultilayerGraphs.jl
index 66f4004..cdac590 100644
--- a/src/MultilayerGraphs.jl
+++ b/src/MultilayerGraphs.jl
@@ -146,7 +146,8 @@ export
     δ_1,
     δ_2,
     δ_3,
-    δ_Ω
+    δ_Ω,
+    havel_hakimi_graph_generator_unsafe
     # tensorfacoriazations.jl
 
 using Base, InteractiveUtils, IterTools, SimpleTraits, Bijections
diff --git a/src/utilities.jl b/src/utilities.jl
index 356ba14..c565da8 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -587,3 +587,47 @@ function havel_hakimi_(empty_graph::SimpleGraph, degree_sequence::Vector{<:Integ
     ##n = length(degree_sequence)
     ##mvs_degree_dict = Dict(mv => deg for (mv, deg) in zip(mvs, degree_sequence) if deg > 0)
 end
+
+
+"""
+    havel_hakimi_(empty_graph::SimpleGraph, degree_sequence::Vector{<: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;
+2. sort the remaining nodes by degree in decreasing order;
+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_graph_generator_unsafe(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.")) =#
+
+    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_ds[v] -= 1
+            println("here")
+        end
+    end
+    
+    
+end

From eecc67e48c496085a2d26fddca66411aa9c0d58a Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Wed, 14 Dec 2022 22:13:26 +0100
Subject: [PATCH 02/17] Create abstractmultilayerugraph_development.jl

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 mains/abstractmultilayerugraph_development.jl | 381 ++++++++++++++++++
 1 file changed, 381 insertions(+)
 create mode 100644 mains/abstractmultilayerugraph_development.jl

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

From 5795a3140a8ff88e441561298d89e5b24dcc6c8d Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Wed, 14 Dec 2022 22:21:46 +0100
Subject: [PATCH 03/17] Implement Havel-Hakimi

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/MultilayerGraphs.jl |  2 +-
 src/utilities.jl        | 54 +++++++++++++++++++++++++++++------------
 2 files changed, 39 insertions(+), 17 deletions(-)

diff --git a/src/MultilayerGraphs.jl b/src/MultilayerGraphs.jl
index cdac590..ffd6994 100644
--- a/src/MultilayerGraphs.jl
+++ b/src/MultilayerGraphs.jl
@@ -147,7 +147,7 @@ export
     δ_2,
     δ_3,
     δ_Ω,
-    havel_hakimi_graph_generator_unsafe
+    havel_hakimi_graph_generator
     # tensorfacoriazations.jl
 
 using Base, InteractiveUtils, IterTools, SimpleTraits, Bijections
diff --git a/src/utilities.jl b/src/utilities.jl
index c565da8..50182e8 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
 """
@@ -590,7 +590,7 @@ end
 
 
 """
-    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;
@@ -601,15 +601,7 @@ Returns a simple graph with a given finite degree sequence of non-negative integ
 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_graph_generator_unsafe(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.")) =#
+function havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer})
 
     graph = SimpleGraph(length(degree_sequence))
     v_ds = OrderedDict(i => deg for (i,deg) in enumerate(degree_sequence))
@@ -618,16 +610,46 @@ function havel_hakimi_graph_generator_unsafe(degree_sequence::Vector{<:Integer})
     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_ds[v] -= 1
-            println("here")
         end
     end
     
+    return graph
     
 end
+
+
+"""
+    kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_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;
+2. sort the remaining nodes by degree in decreasing order;
+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 kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_sequence::AbstractVector{<:Integer})
+
+    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)
+        all(collect(values(v_ds))[1:s] .> 0) || throw(ErrorException("The provided `degree_sequence` is not graphical"))
+        for v in collect(keys(v_ds))[1:s]
+            add_edge!(graph, S, v)
+            v_ds[v] -= 1
+        end
+    end
+    
+    return graph
+    
+end
\ No newline at end of file

From ce24f64046a975de5dda5eb748904dad6acf93b8 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Wed, 14 Dec 2022 22:54:37 +0100
Subject: [PATCH 04/17] Comment the Havel Hakimi graph generator

Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
---
 src/utilities.jl | 38 +++++++++++++++++++++-----------------
 1 file changed, 21 insertions(+), 17 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 50182e8..6d09ef6 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -578,7 +578,7 @@ function havel_hakimi_(empty_graph::SimpleGraph, degree_sequence::Vector{<:Integ
     # 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)."))
+    isgraphical(degree_sequence) || throw(ArgumentError("The degree sequence (degree_sequence) is invalid: it must be graphical (i.e. realizable by 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
@@ -598,30 +598,34 @@ 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)
+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_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))
-    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)
-        all(collect(values(v_ds))[1:s] .> 0) || throw(ErrorException("The provided `degree_sequence` is not graphical"))
-        for v in collect(keys(v_ds))[1:s]
-            add_edge!(graph, S, v)
-            v_ds[v] -= 1
+    # 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 non-negative 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
 
-
 """
     kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_sequence::AbstractVector{<:Integer})
 

From 55a242b25c9c30833f0c973d413681e30edd4a30 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Wed, 14 Dec 2022 23:22:20 +0100
Subject: [PATCH 05/17] Update utilities.jl

Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
---
 src/utilities.jl | 31 -------------------------------
 1 file changed, 31 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 6d09ef6..75edf8e 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -558,37 +558,6 @@ 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})
-
-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;
-2. sort the remaining nodes by degree in decreasing order;
-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 by 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)
-end
-
-
 """
     havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer})
 

From c604cd0a420e80c283b63d6b02f5d2778fe1e0f6 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 00:20:26 +0100
Subject: [PATCH 06/17] Update utilities.jl

Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
---
 src/utilities.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 75edf8e..eeb042c 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -583,7 +583,7 @@ function havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer}
         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 non-negative values
+        # 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]

From 0c85edfb8709b7162e13de99d9c4f6c20374a3cc Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 01:27:10 +0100
Subject: [PATCH 07/17] Implement Kleitman-Wang algorithm (still broken)

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/MultilayerGraphs.jl |  3 +-
 src/utilities.jl        | 74 ++++++++++++++++++++++++++++++++++++-----
 2 files changed, 67 insertions(+), 10 deletions(-)

diff --git a/src/MultilayerGraphs.jl b/src/MultilayerGraphs.jl
index ffd6994..4e6a589 100644
--- a/src/MultilayerGraphs.jl
+++ b/src/MultilayerGraphs.jl
@@ -147,7 +147,8 @@ export
     δ_2,
     δ_3,
     δ_Ω,
-    havel_hakimi_graph_generator
+    havel_hakimi_graph_generator,
+    kleitman_wang_graph_generator
     # tensorfacoriazations.jl
 
 using Base, InteractiveUtils, IterTools, SimpleTraits, Bijections
diff --git a/src/utilities.jl b/src/utilities.jl
index 75edf8e..78a0695 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -595,6 +595,50 @@ function havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer}
     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
+            a < b
+        end
+    end
+
+    return false
+
+end
+
 """
     kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_sequence::AbstractVector{<:Integer})
 
@@ -609,20 +653,32 @@ Returns a simple graph with a given finite degree sequence of non-negative integ
 """
 function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_sequence::AbstractVector{<:Integer})
 
-    graph = SimpleGraph(length(degree_sequence))
-    v_ds = OrderedDict(i => deg for (i,deg) in enumerate(degree_sequence))
+    length(indegree_sequence) == length(outdegree_sequence) || throw(ArgumentError("The provided `indegree_sequence` and `outdegree_sequence` must be of the dame length."))
+
+    graph = SimpleDiGraph(length(indegree_sequence))
+    A_B = zip(deepcopy(indegree_sequence), deepcopy(outdegree_sequence))
+    v_ds = OrderedDict(i => tup for (i,tup) in enumerate(A_B))
     
     
-    while(any(values(v_ds) .!= 0))
-        v_ds = OrderedDict(sort(collect(v_ds), by = last , rev = true))
-        S,s = popfirst!(v_ds)
-        all(collect(values(v_ds))[1:s] .> 0) || throw(ErrorException("The provided `degree_sequence` is not graphical"))
-        for v in collect(keys(v_ds))[1:s]
+    while(any(vcat(values(v_ds)...) .!= 0 ))
+
+        v_ds = OrderedDict(sort(collect(v_ds), by = last, lt = lexicographical_order_ntuple , rev = true))
+
+        @debug "", v_ds
+
+        S, (a_i, b_i) = vcat(popfirst!(v_ds)...)
+
+        # all(vcat(collect(values(v_ds))[1:b_i]...) .> 0) || throw(ErrorException("The degree sequences are not digraphical."))
+
+        @debug "", S, a_i, b_i
+
+        for (v,degs) in collect(v_ds)[1:b_i]
             add_edge!(graph, S, v)
-            v_ds[v] -= 1
+            v_ds[v] = (degs[1]-1, degs[2])
         end
+
     end
     
     return graph
     
-end
\ No newline at end of file
+end

From 2ff4b926c15b4bbffa10183295c3f823c0b15a28 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 18:41:57 +0100
Subject: [PATCH 08/17] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

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"

From 6472e2f6b18f0a4960ed4f3a3bb0da81b0ed542b Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 22:04:05 +0100
Subject: [PATCH 09/17] Implement kleitman-wang (still broken)

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/utilities.jl | 27 ++++++++++++++++++++-------
 1 file changed, 20 insertions(+), 7 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 1b84293..162f109 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -581,7 +581,7 @@ function havel_hakimi_graph_generator(degree_sequence::AbstractVector{<:Integer}
     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 
+        # 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."))
@@ -631,7 +631,7 @@ 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
-            a < b
+            return a < b
         end
     end
 
@@ -660,17 +660,30 @@ function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integ
     v_ds = OrderedDict(i => tup for (i,tup) in enumerate(A_B))
     
     
-    while(any(vcat(values(v_ds)...) .!= 0 ))
+    while(any(Iterators.flatten(values(v_ds)) .!= 0 ))
 
-        v_ds = OrderedDict(sort(collect(v_ds), by = last, lt = lexicographical_order_ntuple , rev = true))
+        
 
-        @debug "", v_ds
+        v_ds = OrderedDict(sort(collect(v_ds), by = last, lt = lexicographical_order_ntuple , rev = true))
 
-        S, (a_i, b_i) = vcat(popfirst!(v_ds)...)
+        @debug "" collect(Iterators.flatten(values(v_ds)))
+
+        # @debug "", v_ds
+        S, (a_i, b_i) = 0,(0,0)
+        for (_S, (_a_i, _b_i)) in collect(deepcopy(v_ds))
+            @debug "" _S _a_i _b_i
+            if _b_i != 0
+                @debug "" _b_i
+                S, a_i, b_i = (_S,_a_i, _b_i) #Iterators.flatten(popfirst!(v_ds)...)
+                delete!(v_ds, _S )
+                break
+            end
+        end
+                
 
         # all(vcat(collect(values(v_ds))[1:b_i]...) .> 0) || throw(ErrorException("The degree sequences are not digraphical."))
 
-        @debug "", S, a_i, b_i
+        # @debug "", S, a_i, b_i
 
         for (v,degs) in collect(v_ds)[1:b_i]
             add_edge!(graph, S, v)

From 6946337935909fe9362cbae129b729ea50a1b1d4 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 22:57:45 +0100
Subject: [PATCH 10/17] Fix kleitman-wang

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/utilities.jl | 64 +++++++++++++++++++++++-------------------------
 1 file changed, 31 insertions(+), 33 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 162f109..b3f6710 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -568,7 +568,7 @@ Returns a simple graph with a given finite degree sequence of non-negative integ
 
 ## 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).
+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
@@ -642,53 +642,51 @@ end
 """
     kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integer},outdegree_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;
-2. sort the remaining nodes by degree in decreasing order;
-3. repeat the procedure.
+Returns a simple directed graph with given finite in-degree and out-degree sequences of non-negative integers generated via the Kleitman-Wang algorithm.
 
 ## 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)
+- [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})
+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))
-    A_B = zip(deepcopy(indegree_sequence), deepcopy(outdegree_sequence))
-    v_ds = OrderedDict(i => tup for (i,tup) in enumerate(A_B))
-    
+    # 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))
     
-    while(any(Iterators.flatten(values(v_ds)) .!= 0 ))
-
-        
-
-        v_ds = OrderedDict(sort(collect(v_ds), by = last, lt = lexicographical_order_ntuple , rev = true))
-
-        @debug "" collect(Iterators.flatten(values(v_ds)))
-
-        # @debug "", v_ds
-        S, (a_i, b_i) = 0,(0,0)
-        for (_S, (_a_i, _b_i)) in collect(deepcopy(v_ds))
-            @debug "" _S _a_i _b_i
+    # 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
-                @debug "" _b_i
-                S, a_i, b_i = (_S,_a_i, _b_i) #Iterators.flatten(popfirst!(v_ds)...)
-                delete!(v_ds, _S )
+                i, a_i, b_i = (_i,_a_i, _b_i)
+                delete!(vertices_degrees_dict, _i)
                 break
             end
         end
                 
+         # Connect the vertex foudn 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
 
-        # all(vcat(collect(values(v_ds))[1:b_i]...) .> 0) || throw(ErrorException("The degree sequences are not digraphical."))
-
-        # @debug "", S, a_i, b_i
+        # 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."))
 
-        for (v,degs) in collect(v_ds)[1:b_i]
-            add_edge!(graph, S, v)
-            v_ds[v] = (degs[1]-1, degs[2])
-        end
+        # Reinsert the vertex, with zero outdegree
+        vertices_degrees_dict[i] = (a_i, 0)
 
     end
     

From 42dae7c859fba9362f4592a64fab22a211f34a00 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Thu, 15 Dec 2022 23:33:15 +0100
Subject: [PATCH 11/17] Add better description to kleitman-wang algorithm

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/utilities.jl | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index b3f6710..15aa4df 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -642,7 +642,11 @@ 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.
+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-oudegree 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-oudegree pairs are of the form (0.0).
 
 ## References
 - [Wikipedia](https://en.wikipedia.org/wiki/Kleitman%E2%80%93Wang_algorithms)

From 8eef50bd2d5fb556d78f2f9a22ba7bbe3c3df0b3 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 13:04:32 +0100
Subject: [PATCH 12/17] Integrate havel-hakimi and kleitman-wang with the
 configuration model-like constructors

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/multilayerdigraph.jl          | 5 ++++-
 src/multilayergraph.jl            | 5 ++++-
 src/utilities.jl                  | 2 +-
 test/abstractmultilayerdigraph.jl | 2 +-
 test/abstractmultilayerugraph.jl  | 2 +-
 5 files changed, 11 insertions(+), 5 deletions(-)

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 15aa4df..6ca7e2a 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -687,7 +687,7 @@ function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integ
         end
 
         # 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."))
+        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)
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)

From 376bd48e91a4dc8beef070b2ca061dded44da4cd Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 15:29:51 +0100
Subject: [PATCH 13/17] Temporarily enable tests on graph_realization branch

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 .github/workflows/CI.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 1534604..2058837 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,7 +3,7 @@ on:
   push:
     branches:
       - main
-      #- dev
+      - graph_realization
     tags: '*'
   pull_request:
 concurrency:

From 1f1f0fa8f0a935e6cfb3eb85b1a0ee72b35db72e Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 15:55:30 +0100
Subject: [PATCH 14/17] Update index.md and README.md

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 README.md         | 1 -
 docs/src/index.md | 1 -
 2 files changed, 2 deletions(-)

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

From 93bd1eaf3226feda196e89724b9402f911a2964f Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 15:57:14 +0100
Subject: [PATCH 15/17] Update utilities.jl

---
 src/utilities.jl | 22 ++++++----------------
 1 file changed, 6 insertions(+), 16 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 6ca7e2a..21c8c3e 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -598,7 +598,7 @@ end
 """
     lexicographical_order_lt(A::Vector{T}, B::Vector{T}) where T
 
-The less than (lt) function that implements lexicographical order
+The less than (lt) function that implements lexicographical order.
 
 See [Wikipedia](https://en.wikipedia.org/wiki/Lexicographic_order).
 """
@@ -611,13 +611,11 @@ function lexicographical_order_lt(A::Union{Vector{T},NTuple{N,T}}, B::Union{Vect
     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
 
 """
@@ -643,10 +641,10 @@ 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-oudegree pairs in lexicographical order;
+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-oudegree pairs are of the form (0.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)
@@ -664,12 +662,10 @@ function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integ
     # 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))
@@ -678,22 +674,16 @@ function kleitman_wang_graph_generator(indegree_sequence::AbstractVector{<:Integ
                 delete!(vertices_degrees_dict, _i)
                 break
             end
-        end
-                
-         # Connect the vertex foudn above to other nodes of highest degree 
+        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
+end
\ No newline at end of file

From 014f409090dda2d3bb5ffbf809e29a44a9c4dcb3 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 16:16:35 +0100
Subject: [PATCH 16/17] Comment old graph realization methods

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 src/utilities.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utilities.jl b/src/utilities.jl
index 21c8c3e..2b93b2c 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -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 =#
 
 
 """

From 3e58f1fd65bdd7e499501907eb3e1b1e5bfa62a7 Mon Sep 17 00:00:00 2001
From: "Interdisciplinary Physics Team (InPhyT)"
 <63546437+InterdisciplinaryPhysicsTeam@users.noreply.github.com>
Date: Fri, 16 Dec 2022 16:18:49 +0100
Subject: [PATCH 17/17] Update CI.yml

Co-Authored-By: Pietro Monticone <38562595+pitmonticone@users.noreply.github.com>
Co-Authored-By: Claudio Moroni <43729990+ClaudMor@users.noreply.github.com>
---
 .github/workflows/CI.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 2058837..5c07a51 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -3,7 +3,6 @@ on:
   push:
     branches:
       - main
-      - graph_realization
     tags: '*'
   pull_request:
 concurrency: