implementation on Hexagonal graph
Closed this issue · 10 comments
I would like to Instantiate an hexagonal graph, so to implement it as a graph type, transforming it into an adjacency matrix (or graph = nx.to_numpy_array(G) if it works better).
Should I use the decorator : @decorator def _nodes_or_number(func_to_be_decorated, *args, **kw) for this matter, and add simply the def hexagonal_lattice_graph(m, n, periodic=False, with_positions=True, create_using=None) from the networkx library?
How about this:
import networkx as nx
from kuramoto import Kuramoto
g = nx.hexagonal_lattice_graph(2, 2)
adj_mat = nx.to_numpy_array(g)
kura = Kuramoto(n_nodes=g.number_of_nodes())
activity = kura.run(adj_mat=adj_mat)Great; so It does yield an order parameter vs coupling constant plot, for an hexagonal grid (as well as a rectangular one which I tried); however, if I take more than 200 nodes for an hexagonal grid, it keeps running endlessly; for the rectangular plot, it looks like the oscillators do not synchronize, even if I take > 2000 nodes (something which for an hexagonal grid seems to be impossible to set).
Is there a way I can fix it to make it run with an hexagonal grid with more than >1000 nodes in an acceptable time interval?
I don't know exactly what you mean by "it runs endlessly". Did you set an order parameter goal to be reached or how does it come that it does not stop? Did you try incrementing the coupling?
It'd be great if you could post a little code snippet so that I can reproduce your issue.
I tried a a bigger range for the coupling constant after seeing that for m, n <= 10 I still couldn't get the oscillators to synchronise
I considered 3 cases: a cyclic, a rectangular and an hexagonal graph, generated using the networkx library.
I chose the natural frequencies for all of them localised at 0 and with "scale=1"; I also chose a quite big range for the coupling, in order to observe the synchronisation of the oscillators at some point (as it takes a stronger coupling depending on the configuration/adjacency matrix).
I also noticed that if I chose a random position of theta as " 2 * np.pi * np.zeros((self.n_nodes,)) " (rather than np.random.random(size=self.n_nodes, for in the kuramoto library, or class Kuramoto), the oscillators synchronise more efficiently; which seems a bit odd as there should not be any difference.
Considering the initial conditions chosen by me, for the three above cites cases I noticed:
-
Cyclic graph, up to 100 nodes: the critical coupling value seems to be bigger than 200, and I am not sure if it makes sense (if I consider the plot of the order parameter vs coupling K)
-
Rectangular and hexagonal grid: up to 10 nodes: the critical coupling constant seems to be bigger than 20; for a number of nodes 100 times 100 (for instance), the synchronisation takes an enormous not realistic amount of time (not excluding I missed something and did not implement it properly).
Also, as for the hexagonal grid (using the networkx library to generate it) we get am 'm x n' rectangular grid of hexagons, which may be what makes the situation even more complicated, but I still fail to understand why even generating an hexagonal lattice with 240 sites (by taking m =10, n = 10), it the code runs without reaching synchronisation within an acceptable time interval (of a few minutes, rather than days) -
I wanted to implement a Runge Kutta method of order 4 method, it looks like the method used is Euler?
I was thinking of using RK45 from the (scipy.integrate)
(in the 'kuramoto' timeseries =odeint(self.derivative, angles_vec, t, args=(adj_mat,))), but I am not sure about how to implement a sensible t_bound (maybe you can suggest the implementation of the integrate.RK45?
the snippet for 1) 2) : (if you could suggest what I am missing /doing wrong):
graph_nx = nx.cycle_graph(n= 100, create_using=None)
graph = nx.to_numpy_array(graph_nx)
kura = Kuramoto(coupling=5, dt=0.01, T=20, n_nodes=len(graph))
act_mat = kura.run(adj_mat=graph)
N = graph_nx.number_of_nodes()
oscillators = 100
Rectangular grid:
RG = nx.grid_2d_graph(m,n)
RGN = nx.convert_node_labels_to_integers(RG)
NR = RGN.number_of_nodes()
pos = {(x,y):(y,-x) for x,y in RG.nodes()}
#pos = nx.spring_layout(RG, seed=42, iterations=100)
adj_mat = nx.to_numpy_array(RG)
kura = Kuramoto(n_nodes=NR)
oscillators = NR
Hexagonal grid:
m = 10
n = 10 #should be possible to generate and have the nodes synchronised without probles as it is a finite grid
#isPBC = True
G = nx.hexagonal_lattice_graph(m,n)
GN = nx.convert_node_labels_to_integers(G)
N = GN.number_of_nodes()
pos = nx.spring_layout(GN, seed=42, iterations=100)
adj_mat = nx.to_numpy_array(G)
kura = Kuramoto(n_nodes=N)
oscillators = N
coupling_vals = np.linspace(0, 50.0, 100)
runs = []
for coupling in coupling_vals:
kura = Kuramoto(coupling=coupling, dt=0.1, T=2500, n_nodes=oscillators)
kura.natfreqs = np.random.normal(loc=0, scale=1, size=oscillators)
act_mat = kura.run(adj_mat=adj_mat)
runs.append(act_mat)
for i, coupling in enumerate(coupling_vals):
r_mean = np.mean([kura.phase_coherence(vec)
for vec in runs_array[i, :, -1000:].T]) # mean over last 1000 steps
plt.plot(coupling, r_mean, 'o')
QUESTION: in one of Strogatz's papers (https://static1.squarespace.com/static/5436e695e4b07f1e91b30155/t/544525a8e4b0b8e2e8230fa3/1413817768995/from-kuramoto-to-crawford.pdf) the critical coupling is taken as 2/(pi * g(0)) which is g at omega = 0, fixing for instance omega to be a Gaussian; how can I relate to the predicted critical coupling defined in your code?
As I pointed out here (https://github.com/fabridamicelli/kuramoto#kuramoto-model-201), bear in mind the following points:
- synchronization is not guaranteed, ie might not occur at all.
- Partial synchronization is a possible outcome.
- The higher the dimension of the lattice on which the oscillators are embedded, the easier it is to synchronize. For example, there isn't any good synchronization in one dimension, even with strong coupling. In two dimensions it is not clear yet. From 3 dimensions on, the model starts behaving more like the mean field prediction.
If I understand correctly, the graphs that you're dealing with are all embedded in two dimensions, so that might be an important factor.
Critical coupling value in my example was calculated after this https://iopscience.iop.org/article/10.1088/0143-0807/29/1/015.
Regarding the integration methods, I'm afraid I don't know them enough to give you advice.
But before trying that, I would rather go for trying to break the two dimensionality of your graphs, for example, adding some random links, and see if that allows them to synchronize better
- "synchronization is not guaranteed, ie might not occur at all." in the case of the hexagonal lattice, for a number of nodes above 200, it does converge if we wait long enough, but that long enough is extremely long compared to the performance I'd like to reach, which however, may not be doable within a shorter time.
What about choosing some convergence as stopping after some value is achieved (i.e. implemeting a while loop as "while (error > stop or error step < averaging) and i < 1000 etc..." for instance), so stopping when it is smaller than a fixed error (for instance 10^-10) and letting the coupling loop starting from the biggest value to the smallest one (e.g. from 50 down to 0; I'm trying the latter, with little success in shortening the time)?
Any suggestion as to how I could implement that?
-as you pointed out, I'm interested in a 2D lattice; I am also wondering if instead of the adjacency matrix one uses the Incidence matrix it'd help in terms of computation time? (I would need some hints as to how to change that in the Kuramoto class)
-
""Partial synchronization is a possible outcome"": is that what I observe for an hexagonal lattice (for ex.), when choosing the initial angles as defined by you (i.e 2 * np.pi * np.random.random(size=self.n_nodes)?
I set them to np.zeros((self.n_nodes,)), and it seems to fully synchronise in most cases, which makes me wonder why, as there shouldn't be any difference? -
Also, I was thinking of storing the last values after synchronisation to use them as the initial ones, would you have any suggestions as to how to implement it?
[thanks for the link, I also saw that paper, but I used it for constructing an all-to-all model; I am also working on another code but without classes and objects, more in the functional way, so that I can have a comparison].
" I would rather go for trying to break the two dimensionality of your graphs, for example, adding some random links, and see if that allows them to synchronise better" :
Does the periodicity condition count for that too (that is, if it is periodic would the "additional" edges count as additional links?)
When the synchronization happens (if it happens at all) it is typically an abrubt transition. So I am not sure how fruitful it is to use a threshold. But yes, you could implement a loop with a stop criterion based on the synchronization level and some "patience" number of steps to avoid runnning into an infinite loop.
I don't think that changing the representation of the connectivity would help. The model operates at the node level, which one can read directly from rows/columns of the adjacency matrix.
I am not sure what you mean with periodicity. I just meant simply adding links between nodes taken at random. In that way one allows to "jump" over the grid neighbors, which, again, breaks the two-dimensionality of your graph.
"ut yes, you could implement a loop with a stop criterion based on the synchronization level and some "patience" number of steps to avoid runnning into an infinite loop." aisn't this patience number of steps (and error) for the convergence criterium given by the (optional) parameters "atol, rtol" of the solver "scipy.integrate.odeint" ? (ref. https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html)
"I don't think that changing the representation of the connectivity would help. The model operates at the node level, which one can read directly from rows/columns of the adjacency matrix." Also from a simpler code I wrote I noticed the computation time was very long even with a different representation of the connectivity.
"I am not sure what you mean with periodicity." that is the option in networkx to make a periodic grid by joining the boundary vertices. Not sure how that could help.
I also tried to add the Runge Kutta 4th order method to solve the Kuramoto equation but (implementation aside which may be arguably right) it does not help getting at synchronisation; according to other sources (where Fortran or Matlab was used) there should be synchrony for an hexagonal grid. I am not sure what the issue is in my case.