Creating a layered graph (aka messing with property maps)

Hi all,
I'm trying to analyze a graph for which I have multiple layers and I'm finding difficulties in creating the graph. I have two graphs with the same vertex for which I define a edge property (the layer) like this:

(g_atac, dist_atac) = gt.generate_knn(A_data, k=n_neighbors)
(g_rna, dist_rna) = gt.generate_knn(R_data, k=n_neighbors)

layer_atac = g_atac.new_edge_property('int')
for e in g_atac.edges():
    layer_atac[e] = 1
layer_rna = g_rna.new_edge_property('int')
for e in g_rna.edges():
    layer_rna[e] = -1
g_atac.edge_properties['layer'] = layer_atac
g_rna.edge_properties['layer'] = layer_rna

I then merge the graphs by graph union, first defining a node mapping:

rna_mappings = g_rna.new_vertex_property("int")
for x in range(g_atac.num_vertices()):
    rna_mappings[x] = x

and then performing the actual union:

gu = gt.graph_union(g_atac, g_rna, intersection=rna_mappings, internal_props=True)

This indeed creates the final graph, but I noticed that instead of two layers I have three:

np.unique(gu.edge_properties['layer'].a)
PropertyArray([-1, 0, 1], dtype=int32)

Looking back at the start graphs I have zeros in the edge_properties as well

np.unique(g_rna.edge_properties['layer'].a)
PropertyArray([-1, 0], dtype=int32)

np.unique(g_atac.edge_properties['layer'].a)
PropertyArray([0, 1], dtype=int32)

and, similarly, the elements in the edge property vector are much higher than the number of edges:

print((len(g_rna.ep['layer'].a), g_rna.num_edges()))
(156672, 14317)

which is different from what I can get from this:

g_rna.get_edges([g_rna.ep['layer']]).shape
(14317, 3)

I'm evidently messing with (internal) properties and I'm clueless, any advice?

As a side note, if I try to create a LayeredBlockState from any of the initial graphs I only have one block:

gt.NestedBlockState(g_atac, base_type=gt.LayeredBlockState, state_args=dict(ec=g_atac.ep.layer, layered=True))
<NestedBlockState object, with base <LayeredBlockState object with 1 blocks, 1 edge covariates, degree-corrected, for graph <Graph object, undirected, with 1088 vertices and 16529 edges, 1 internal edge property, at 0x7ffb4eaa17c0>, at 0x7ffb4eaf0610>, and 12 levels of sizes [(1088, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1)] at 0x7ffb4eaf0340>

I noticed the whole thing because I have instead 3 layers in the union graph:

gt.NestedBlockState(gu, base_type=gt.LayeredBlockState, state_args=dict(ec=gu.ep.layer, layered=True))
<NestedBlockState object, with base <LayeredBlockState object with 1 blocks, 3 edge covariates, degree-corrected, for graph <Graph object, undirected, with 1088 vertices and 30846 edges, 1 internal edge property, at 0x7ffb4ea9dac0>, at 0x7ffb4eaf0730>, and 12 levels of sizes [(1088, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1), (1, 1)] at 0x7ffb4eaf0550>

and some of the edges with layer 0 in the union exist in the original graphs with assigned layer 1 or -1

I may have found a way to go, but I'd like to know what I was doing wrong before.
Now I don't create a new graph by union, but this (as vertex numbers and indexes are the same):

gu = gt.Graph(directed=False)
gu.add_vertex(g_atac.num_vertices())
layer = gu.new_edge_property('int')
for e in g_rna.edges():
    ne = gu.add_edge(e.source(), e.target())
    layer[ne] = -1
for e in g_atac.edges():
    ne = gu.add_edge(e.source(), e.target())
    layer[ne] = 1
gu.ep['layer'] = layer

this works and the final graph has 2 layers and no extra zeros

As is explained in the documentation, edge indexes need not to be
contiguous (differently from vertex indexes, which are always contiguous).

When accessing the property map values via the ".a" property, this
returns an array indexed by the edge index, which will contain values
for non-existing edges if the indexes are not contiguous. This is what
you are seeing.

If you want to see only the property values for existing edges, you
should get a filtered array with the ".fa" property.

Best,
Tiago