Kind of pooling to use with node features representing longitude/lattitude

[StefanBloemheuvel]

Hi creators of Spektral,

First of all, thank you for creating this amazing package.
I am working on combining graph neural networks with time series data. In here, I have as node features for the graph convolutional layer: longitude and lattitude of each node in the network. As the graph input i have the normalized adjacency matrices. The other inputs are my time series data for each node which go in the cnn.

What would be a correct way to use pooling in this scenario? The decision confuses me since it does not make sense to average, sum or max the node features, since a higher or lower longitude or lattitude is not better than a lower longitude lattitude). Also, the relevance of a node is only made clear after the graph convolutional layers where i append a classical cnn to process the time series data. So using attention layers for pooling makes no sense to me as well, since you need the time series to know which node has some interesting things going on.

[danielegrattarola]

Hi,

Depending on how you structured your network (feel free to post a model.summary() of your architecture) you should not be pooling the node features without having processed them with the GNN first.
Global pooling is similar to the flattening operation in classical CNNs, so you only should do it when you're ready to feed your "flattened" representation to a MLP or something along those lines.

If you do process the features with a GNN, then the overall network will learn to compute features that make sense when summed (or averaged, or maxed) together. It's not a matter of being "better", it's a matter of being informative enough to solve your task.

I am not sure I understand how you process your time series, so I'd need more details to give a more precise answer, but the idea is more or less the same: first you run your convolutions, then pool, and finally run your classification/regression.

On a more practical note: the choice of global pooling may depend on your task, but sum pooling is usually a safe bet.

Cheers

[StefanBloemheuvel]

I have tried a lot of different configurations and some things are starting to notice me. I have a kinda unique dataset of earthquake measuring stations and they each measure several earthquakes (time series data). I have created a weighed graph of all the pairwise distances between the stations. However, not each earthquake is measured by all the stations. So some earthquakes are measured by 60+ stations and others only by 30+/-.

Therefore, in order to find a nice balance, I had to select for each earthquake a random sample of 50 stations to use in the cnn (since they always need the same shape as input) and also for the graph. This means that if an earthquake was measured by less than 50 stations, the missing number of stations will be filled in by 'fake' stations that only have time series readings of zeros (and no stations locations as node features for the graph).

I only found out recently that having such 'fake' stations really messes with the results. This is due to an error in my code that allowed some edges to still be created between the 'real' stations and the 'fake' stations. So i thought the 'fake' stations had no connections in the graph, but they did. And when I thought I would get better results when the faulty edges were removed, my results were 10x worse in predicting some features of the earthquakes. So it seems that having these connections to the fake nodes is necessary.

E.g., a generated example of a weighted adjacency matrix of me would like this:

import numpy as np
import networkx as nx

G = nx.Graph()

G.add_edge("a", "b", weight=0.6)
G.add_edge("a", "c", weight=0.2)
G.add_edge("c", "d", weight=0.2)
G.add_edge("c", "e", weight=0.6)
G.add_edge("c", "f", weight=0.6)
G.add_edge("a", "d", weight=0.3)

G.add_node("g") # <-- this would be a station without any connections

pos = nx.spring_layout(G, seed=7)
weights_calculated = nx.get_edge_attributes(G,'weight').values()
nx.draw(G, pos, 
        width=list(weights_calculated),
        with_labels=True)

ax = plt.gca()
ax.margins(0.08)
plt.axis("off")
plt.tight_layout()
plt.show()

adj = nx.adjacency_matrix(G).todense()
print(adj)

[[0.  0.6 0.2 0.3 0.  0.  0. ]
 [0.6 0.  0.  0.  0.  0.  0. ]
 [0.2 0.  0.  0.2 0.6 0.6 0. ]
 [0.3 0.  0.2 0.  0.  0.  0. ]
 [0.  0.  0.6 0.  0.  0.  0. ]
 [0.  0.  0.6 0.  0.  0.  0. ]
 [0.  0.  0.  0.  0.  0.  0. ]]

In order to make this adjacency matrix applicable for a graph neural network, we for example would apply the gcn_filter to it, which returns this weighted matrix:

adj_new = gcn_filter(adj, symmetric = True)
adj_new = np.round(adj_new, 2)
adj_new
array([[0.48, 0.33, 0.09, 0.17, 0.  , 0.  , 0.  ],
       [0.33, 0.62, 0.  , 0.  , 0.  , 0.  , 0.  ],
       [0.09, 0.  , 0.38, 0.1 , 0.29, 0.29, 0.  ],
       [0.17, 0.  , 0.1 , 0.67, 0.  , 0.  , 0.  ],
       [0.  , 0.  , 0.29, 0.  , 0.62, 0.  , 0.  ],
       [0.  , 0.  , 0.29, 0.  , 0.  , 0.62, 0.  ],
       [0.  , 0.  , 0.  , 0.  , 0.  , 0.  , 1.  ]])

Notice the 1 in the bottom-right corner still being there, having no other connections to other nodes? I think this messes with the results. Do you perhaps have an explanation for this? And should I then just create a graph with connections to these 'fake' stations with locations set to (0,0), as the only way to model these graphs?

Thanks in advance! :)

[danielegrattarola]

Graph neural networks are designed to take inputs of different sizes, so I would first try that and see if the performance changes.

I would suggest that you convert your data to the spektral.data API so that you can create a custom dataset and use a DisjointLoader to create batches of differently sized graphs.

Alternatively, you can try to set the extra 1s that you have on the diagonal to 0, so that any GNN computation will ignore the station entirely.

I think the way in which you create a graph from the data is extremely specific to the problem, I cannot say a prior if adding extra stations will significantly interfere with the results. On the surface, it doesn't look like it should disrupt the GNN entirely, at most you would expect a small-ish drop in performance.

Cheers

[StefanBloemheuvel]

Thank you for the fast reply! I will first try the second option right away and put the 1's of the 'fake' stations to 0, since that will be most easily tried.

However, the first option you mention sounds most ideal. Will the DisjointLoader solve the problem that the cnn expects the same input shape each time? The data shape for my cnn input is (None, N_stations, 2048, N_channels), so in each earthquake the N_stations will then be different, which I thought was not possible in cnn's?

[danielegrattarola]

Yes, that's one of the cool things of GNNs (also note... the CNN could, in principle, work for differently-sized inputs but the problem is that we need to pack images in a tensor where dimensions are fixed).
In disjoint mode you will simple store your data as a Dataset of Graph objects, and the disjoint loader will take care of the rest.

Have a look at the examples to see how that works.

Cheers

[ecastillot]

ah that clears thing up really wel considering the pooling being similar to flattening, and also the different pooling techniques!

One question though is your formulation on how to process the time series, since I am processing the time series in a seperate cnn and concatenate the features of the cnn with the features of the gcnn before the dense layers.

So my model now in principle is:

Graph layer (input is graph and node features of the sensors)
Graph layer
Pooling

CNN layer (input are the time series of the sensors)
CNN layer
CNN layer
flatten

Concat (pooling, flatten)
Dense
Dense output

So the gcnn does not process the time series here, or is that actually a better method?

Dear @StefanBloemheuvel ,
Could you be so kind as to share your architecture model. It could be really useful for me. I have similar type of data focused in seismology subject.
I want to use CNN layers for time series signal and after concatenate them with location features. Your scheme would help me to make a lot of progress.

Thank you so much.

[StefanBloemheuvel]

Hi @danielegrattarola, I got pretty far with adjusting your example dataset generating .py file to my kind of usecase. It is not a regression task, but now it is an toy usecase inspired from the 'multiple parallel series' paragraph from here. The class generates a set of n nodes, gives it x and y locations by giving the networkx graph a layout, calculates the distance between each node as an edge weight, then generates similar time series as from the website example named above. In that example, they create multiple parallel series and want to predict the next value for each series:

10, 15, 25, .., n
20, 25, 45, .., n
30, 35, 65, .., n

returns

40, 45, 85, .., n

So the graph has the adjacency matrix, the xy locations as node features, the time series as data, and the node labels as y are the next values in the time series.

import matplotlib.pyplot as plt
from spektral.layers import GlobalSumPool
from scipy.spatial.distance import euclidean 
import random
import networkx as nx
import numpy as np 
np.set_printoptions(suppress=True)
from spektral.data import Dataset, DisjointLoader, Graph

class MyGraphMaker(Dataset):
    '''
    We try to predict the next value in a timeseries of length 3 for n nodes.
    data shape for each graph is (forecast_samples, number_of_time_steps, num_of_nodes)
    When length is 9 and n_steps is 3, forecast_samples will be 6.
    So each graph will have (6,3,num_of_nodes)

    The shape of y is (6,num_of_nodes)

    for each node, the cnn should predict the next value in the steps taken. 
    So for node 1 with n_steps = 3 the sequence is e.g., [0,1,2], where the y value is 3. 
    '''

    def __init__(self, n_samples, n_min=5, n_max=30, p = 0.7, n_steps = 3, length=9, **kwargs):
        self.n_samples = n_samples
        self.n_min = n_min
        self.n_max = n_max
        self.p = p
        self.n_steps = n_steps
        self.length = length
        super().__init__(**kwargs)

    def split_sequences(self, sequences, n_steps):
        X, y = list(), list()
        for i in range(len(sequences)):
            # find the end of this pattern
            end_ix = i + n_steps
            # check if we are beyond the dataset
            if end_ix > len(sequences)-1:
                break
            # gather input and output parts of the pattern
            seq_x, seq_y = sequences[i:end_ix, :], sequences[end_ix, :]
            X.append(seq_x)
            y.append(seq_y)
        return np.array(X), np.array(y)


    def read(self):
        def make_graph():

            #determine the number of nodes in the graph
            num_nodes = np.random.randint(self.n_min, self.n_max)

            # create the graph, where the weights in the adjacency matrix represent
            # the euclidean distance between the x,y locations of each node, reverse weighted by 1/distance to map them around [0,1]
            G = nx.connected_watts_strogatz_graph(n = num_nodes, k = 4, p=self.p, tries = 3)
            pos = nx.spring_layout(G,)# you can add a seed here
            
            edge_updates = {}
            for (a,b,_) in G.edges(data=True): 
                distance = euclidean(pos[a], pos[b])
                edge_updates[a,b] = {'weight': 1/distance}
            nx.set_edge_attributes(G, edge_updates)
            adj = np.round(nx.adjacency_matrix(G, weight='weight', nodelist=nx.nodes(G)),2)
            
            # create the time series of length = self.length, where each node has values going from low to high
            # the cnn has to predict the next value in the sequence
            # the starting value is based on the y value of the x,y location, mimicking that that 'matters'.
            timeseries = np.zeros((self.length,num_nodes))
            for i,j in enumerate(pos.items()):
                value = np.round(j[1][1],2)
                timeseriesnode = np.zeros(self.length)
                timeseriesnode = []
                for k in range(0,self.length): 
                    timeseriesnode.append(value)
                    value += 1
                timeseries[:,i] = timeseriesnode
            timeseries, outcome = self.split_sequences(timeseries, self.n_steps)

            # create the node features, which are the x,y locations of the nodes
            positions = np.zeros((num_nodes, 2))
            for q,w in pos.items():
                positions[q,:] = w
            
            return Graph(x=positions, a=adj, data=timeseries, y=outcome)

        # We must return a list of Graph objects
        return [make_graph() for _ in range(self.n_samples)]

dataset = MyGraphMaker(32)#, transforms=NormalizeAdj())
#%%
print(dataset[0])
print(dataset[0], '\n')

F = dataset.n_node_features  # Dimension of node features
n_out = dataset.n_labels  # Dimension of the target
print(F)
print(n_out, '\n')

try:
    print(dataset.data)
except:
    print('not found data', '\n')

print('however, the data can be found when indexing the dataset class:')
print(dataset[0].data.shape)
print(dataset[0].y.shape)

However, I have some trouble in accessing the time series data that I give the MyGraphMaker class by returning:

return Graph(x=positions, a=adj, data=timeseries, y=outcome)

When I try to print dataset.data, it will give an attribute error saying the MyGraphMaker object has no attribute 'data'. How could I tell the class that I would like to retrieve this data variable for each of the graphs in the class? Similar to doing dataset.n_labels?
It is however, possible to retrieve the data by indexing, so dataset[0].data works, but not dataset.data.

Kind regards! :)

Stefan

[danielegrattarola]

That is correct! You can subclass DisjointLoader and re-implement the collate() method.

Let me know if you need help in the process

[StefanBloemheuvel]

Okay I have been digging a bit into the code, but I stumble upon the following error when creating a batch of graphs in DisjointMode (no altering yet).

print(dataset[0])
print(dataset[1])


>>> Graph(n_nodes=27, n_node_features=2, n_edge_features=None, n_labels=27)
>>> Graph(n_nodes=26, n_node_features=2, n_edge_features=None, n_labels=26)

loader_tr = DisjointLoader(dataset_tr, batch_size=batch_size, epochs=epochs, node_level=True)

batch = loader_tr.__next__()
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~/Documents/GitHub/Earthquake_GNN/pythonfiles/disjoint.py in <module>
     25 loader_tr = DisjointLoader(dataset_tr, batch_size=batch_size, epochs=epochs, node_level=True)
     26 
---> 27 batch = loader_tr.__next__()
     28 
     29 

~/anaconda3/envs/earthquake/lib/python3.7/site-packages/spektral_stefan/data/loaders.py in __next__(self)
     81     def __next__(self):
     82         nxt = self._generator.__next__()
---> 83         return self.collate(nxt)
     84 
     85     def generator(self):

~/anaconda3/envs/earthquake/lib/python3.7/site-packages/spektral_stefan/data/loaders.py in collate(self, batch)
    261             # print(np.vstack(y))
    262             # print(y)
--> 263             y = np.vstack(y) if self.node_level else np.array(y)
    264         output = to_disjoint(**packed)
    265 

<__array_function__ internals> in vstack(*args, **kwargs)

~/anaconda3/envs/earthquake/lib/python3.7/site-packages/numpy/core/shape_base.py in vstack(tup)
    281     if not isinstance(arrs, list):
    282         arrs = [arrs]
--> 283     return _nx.concatenate(arrs, 0)
    284 
    285 

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: all the input array dimensions for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 27 and the array at index 1 has size 26

This 27 and 26 come from the graphs having 27 and 26 nodes, and the values in n_labels represent the next value in the time series that each node holds. So having a graph with 210 nodes will mean the n_labels will have a length of 210.

It seems to me that the package only allows a flexible number of nodes, but not labels? But that is strange to me, since graphs having different number of nodes should have different number of labels (in node_level=True). Could perhaps a workaround be that I do not supply the graph dataset the node labels, but try to also pass them through with my still-to-be-made altered DisjointLoader? The n_labels aren't even node labels anyway, they actually are just the to be predicted values from the sequence time series. That is totally different from e.g., predicting the label of the node which is only one value of course.

Kind regards,

Stefan

[danielegrattarola]

DisjointLoader assumes that the labels are of shape (n_nodes, n_labels) and tries to stack them together. The error is due to the fact that you have 1-dimensional labels but the dimension is implicit.

Can you try to convert your labels to be of shape (n_nodes, 1)? It should fix your issues

Cheers

[StefanBloemheuvel]

Hi, I was busy with some other work but back on track with this project again! Yes, indeed reshaping the labels worked, now I will start creating my own Loader. Is it true that I could use this Q&A https://github.com/danielegrattarola/spektral/discussions/141 as a starting-ground for my own usecase?

[danielegrattarola]

Yes, that's definitely a good starting point!