sync.py
1 """!
2 
3 @brief Neural Network: Oscillatory Neural Network based on Kuramoto model
4 @details Implementation based on paper @cite article::syncnet::1, @cite article::nnet::sync::1, @cite inproceedings::net::sync::1.
5 
6 @authors Andrei Novikov (pyclustering@yandex.ru)
7 @date 2014-2018
8 @copyright GNU Public License
9 
10 @cond GNU_PUBLIC_LICENSE
11  PyClustering is free software: you can redistribute it and/or modify
12  it under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  PyClustering is distributed in the hope that it will be useful,
17  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  GNU General Public License for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with this program. If not, see <http://www.gnu.org/licenses/>.
23 @endcond
24 
25 """
26 
27 import math
28 import numpy
29 import random
30 import warnings
31 
32 try:
33  import matplotlib.pyplot as plt
34  import matplotlib.animation as animation
35 except Exception as error_instance:
36  warnings.warn("Impossible to import matplotlib (please, install 'matplotlib'), pyclustering's visualization "
37  "functionality is not available (details: '%s')." % str(error_instance))
38 
39 import pyclustering.core.sync_wrapper as wrapper
40 
41 from pyclustering.core.wrapper import ccore_library
42 
43 from scipy.integrate import odeint
44 
45 from pyclustering.nnet import network, conn_represent, conn_type, initial_type, solve_type
46 from pyclustering.utils import pi, draw_dynamics, draw_dynamics_set, set_ax_param
47 
48 
50  """!
51  @brief Provides services to calculate order parameter and local order parameter that are used
52  for synchronization level estimation.
53 
54  """
55  @staticmethod
56  def calculate_sync_order(oscillator_phases):
57  """!
58  @brief Calculates level of global synchronization (order parameter) for input phases.
59  @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when
60  desynchronization is observed in the network.
61 
62  @param[in] oscillator_phases (list): List of oscillator phases that are used for level of global synchronization.
63 
64  @return (double) Level of global synchronization (order parameter).
65 
66  @see calculate_order_parameter()
67 
68  """
69 
70  exp_amount = 0.0;
71  average_phase = 0.0;
72 
73  for phase in oscillator_phases:
74  exp_amount += math.expm1( abs(1j * phase) );
75  average_phase += phase;
76 
77  exp_amount /= len(oscillator_phases);
78  average_phase = math.expm1( abs(1j * (average_phase / len(oscillator_phases))) );
79 
80  return abs(average_phase) / abs(exp_amount);
81 
82 
83  @staticmethod
84  def calculate_local_sync_order(oscillator_phases, oscillatory_network):
85  """!
86  @brief Calculates level of local synchorization (local order parameter) for input phases for the specified network.
87  @details This parameter is tend 1.0 when the oscillatory network close to local synchronization and it tend to 0.0 when
88  desynchronization is observed in the network.
89 
90  @param[in] oscillator_phases (list): List of oscillator phases that are used for level of local (partial) synchronization.
91  @param[in] oscillatory_network (sync): Instance of oscillatory network whose connections are required for calculation.
92 
93  @return (double) Level of local synchronization (local order parameter).
94 
95  """
96 
97  exp_amount = 0.0;
98  num_neigh = 0.0;
99 
100  for i in range(0, len(oscillatory_network), 1):
101  for j in range(0, len(oscillatory_network), 1):
102  if (oscillatory_network.has_connection(i, j) == True):
103  exp_amount += math.exp(-abs(oscillator_phases[j] - oscillator_phases[i]));
104  num_neigh += 1.0;
105 
106  if (num_neigh == 0):
107  num_neigh = 1.0;
108 
109  return exp_amount / num_neigh;
110 
111 
112 
114  """!
115  @brief Represents output dynamic of Sync.
116 
117  """
118 
119  @property
120  def output(self):
121  """!
122  @brief (list) Returns output dynamic of the Sync network (phase coordinates of each oscillator in the network) during simulation.
123 
124  """
125  if ( (self._ccore_sync_dynamic_pointer is not None) and ( (self._dynamic is None) or (len(self._dynamic) == 0) ) ):
126  self._dynamic = wrapper.sync_dynamic_get_output(self._ccore_sync_dynamic_pointer);
127 
128  return self._dynamic;
129 
130 
131  @property
132  def time(self):
133  """!
134  @brief (list) Returns sampling times when dynamic is measured during simulation.
135 
136  """
137  if ( (self._ccore_sync_dynamic_pointer is not None) and ( (self._time is None) or (len(self._time) == 0) ) ):
138  self._time = wrapper.sync_dynamic_get_time(self._ccore_sync_dynamic_pointer);
139 
140  return self._time;
141 
142 
143  def __init__(self, phase, time, ccore = None):
144  """!
145  @brief Constructor of Sync dynamic.
146 
147  @param[in] phase (list): Dynamic of oscillators on each step of simulation. If ccore pointer is specified than it can be ignored.
148  @param[in] time (list): Simulation time.
149  @param[in] ccore (ctypes.pointer): Pointer to CCORE sync_dynamic instance in memory.
150 
151  """
152 
153  self._dynamic = phase;
154  self._time = time;
155  self._ccore_sync_dynamic_pointer = ccore;
156 
157 
158  def __del__(self):
159  """!
160  @brief Default destructor of Sync dynamic.
161 
162  """
163  if (self._ccore_sync_dynamic_pointer is not None):
164  wrapper.sync_dynamic_destroy(self._ccore_sync_dynamic_pointer);
165 
166 
167  def __len__(self):
168  """!
169  @brief Returns number of simulation steps that are stored in dynamic.
170  @return (uint) Number of simulation steps that are stored in dynamic.
171 
172  """
173  if (self._ccore_sync_dynamic_pointer is not None):
174  return wrapper.sync_dynamic_get_size(self._ccore_sync_dynamic_pointer);
175 
176  return len(self._dynamic);
177 
178 
179  def __getitem__(self, index):
180  """!
181  @brief Indexing of the dynamic.
182 
183  """
184  if (index is 0):
185  return self.time;
186 
187  elif (index is 1):
188  return self.output;
189 
190  else:
191  raise NameError('Out of range ' + index + ': only indexes 0 and 1 are supported.');
192 
193 
194  def allocate_sync_ensembles(self, tolerance = 0.01, indexes = None, iteration = None):
195  """!
196  @brief Allocate clusters in line with ensembles of synchronous oscillators where each synchronous ensemble corresponds to only one cluster.
197 
198  @param[in] tolerance (double): Maximum error for allocation of synchronous ensemble oscillators.
199  @param[in] indexes (list): List of real object indexes and it should be equal to amount of oscillators (in case of 'None' - indexes are in range [0; amount_oscillators]).
200  @param[in] iteration (uint): Iteration of simulation that should be used for allocation.
201 
202  @return (list) Grours (lists) of indexes of synchronous oscillators.
203  For example [ [index_osc1, index_osc3], [index_osc2], [index_osc4, index_osc5] ].
204 
205  """
206 
207  if (self._ccore_sync_dynamic_pointer is not None):
208  ensembles = wrapper.sync_dynamic_allocate_sync_ensembles(self._ccore_sync_dynamic_pointer, tolerance, iteration);
209 
210  if (indexes is not None):
211  for ensemble in ensembles:
212  for index in range(len(ensemble)):
213  ensemble[index] = indexes[ ensemble[index] ];
214 
215  return ensembles;
216 
217  if ( (self._dynamic is None) or (len(self._dynamic) == 0) ):
218  return [];
219 
220  number_oscillators = len(self._dynamic[0]);
221  last_state = None;
222 
223  if (iteration is None):
224  last_state = self._dynamic[len(self._dynamic) - 1];
225  else:
226  last_state = self._dynamic[iteration];
227 
228  clusters = [];
229  if (number_oscillators > 0):
230  clusters.append([0]);
231 
232  for i in range(1, number_oscillators, 1):
233  cluster_allocated = False;
234  for cluster in clusters:
235  for neuron_index in cluster:
236  last_state_shifted = abs(last_state[i] - 2 * pi);
237 
238  if ( ( (last_state[i] < (last_state[neuron_index] + tolerance)) and (last_state[i] > (last_state[neuron_index] - tolerance)) ) or
239  ( (last_state_shifted < (last_state[neuron_index] + tolerance)) and (last_state_shifted > (last_state[neuron_index] - tolerance)) ) ):
240  cluster_allocated = True;
241 
242  real_index = i;
243  if (indexes is not None):
244  real_index = indexes[i];
245 
246  cluster.append(real_index);
247  break;
248 
249  if (cluster_allocated == True):
250  break;
251 
252  if (cluster_allocated == False):
253  clusters.append([i]);
254 
255  return clusters;
256 
257 
258  def allocate_phase_matrix(self, grid_width = None, grid_height = None, iteration = None):
259  """!
260  @brief Returns 2D matrix of phase values of oscillators at the specified iteration of simulation.
261  @details User should ensure correct matrix sizes in line with following expression grid_width x grid_height that should be equal to
262  amount of oscillators otherwise exception is thrown. If grid_width or grid_height are not specified than phase matrix size
263  will by calculated automatically by square root.
264 
265  @param[in] grid_width (uint): Width of the allocated matrix.
266  @param[in] grid_height (uint): Height of the allocated matrix.
267  @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
268  If iternation number is not specified, the last step of simulation is used for the matrix allocation.
269 
270  @return (list) Phase value matrix of oscillators with size [number_oscillators x number_oscillators].
271 
272  """
273 
274  output_dynamic = self.output;
275 
276  if ( (output_dynamic is None) or (len(output_dynamic) == 0) ):
277  return [];
278 
279  current_dynamic = output_dynamic[len(output_dynamic) - 1];
280  if (iteration is not None):
281  current_dynamic = output_dynamic[iteration];
282 
283  width_matrix = grid_width;
284  height_matrix = grid_height;
285  number_oscillators = len(current_dynamic);
286  if ( (width_matrix is None) or (height_matrix is None) ):
287  width_matrix = int(math.ceil(math.sqrt(number_oscillators)));
288  height_matrix = width_matrix;
289 
290  if (number_oscillators != width_matrix * height_matrix):
291  raise NameError("Impossible to allocate phase matrix with specified sizes, amout of neurons should be equal to grid_width * grid_height.");
292 
293 
294  phase_matrix = [ [ 0.0 for i in range(width_matrix) ] for j in range(height_matrix) ];
295  for i in range(height_matrix):
296  for j in range(width_matrix):
297  phase_matrix[i][j] = current_dynamic[j + i * width_matrix];
298 
299  return phase_matrix;
300 
301 
302  def allocate_correlation_matrix(self, iteration = None):
303  """!
304  @brief Allocate correlation matrix between oscillators at the specified step of simulation.
305 
306  @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
307  If iternation number is not specified, the last step of simulation is used for the matrix allocation.
308 
309  @return (list) Correlation matrix between oscillators with size [number_oscillators x number_oscillators].
310 
311  """
312 
313  if (self._ccore_sync_dynamic_pointer is not None):
314  return wrapper.sync_dynamic_allocate_correlation_matrix(self._ccore_sync_dynamic_pointer, iteration);
315 
316  if ( (self._dynamic is None) or (len(self._dynamic) == 0) ):
317  return [];
318 
319  dynamic = self._dynamic;
320  current_dynamic = dynamic[len(dynamic) - 1];
321 
322  if (iteration is not None):
323  current_dynamic = dynamic[iteration];
324 
325  number_oscillators = len(dynamic[0]);
326  affinity_matrix = [ [ 0.0 for i in range(number_oscillators) ] for j in range(number_oscillators) ];
327 
328  for i in range(number_oscillators):
329  for j in range(number_oscillators):
330  phase1 = current_dynamic[i];
331  phase2 = current_dynamic[j];
332 
333  affinity_matrix[i][j] = abs(math.sin(phase1 - phase2));
334 
335  return affinity_matrix;
336 
337 
338  def calculate_order_parameter(self, start_iteration = None, stop_iteration = None):
339  """!
340  @brief Calculates level of global synchorization (order parameter).
341  @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when
342  desynchronization is observed in the network. Order parameter is calculated using following equation:
343 
344  \f[
345  r_{c}=\frac{1}{Ne^{i\varphi }}\sum_{j=0}^{N}e^{i\theta_{j}};
346  \f]
347 
348  where \f$\varphi\f$ is a average phase coordinate in the network, \f$N\f$ is an amount of oscillators in the network.
349 
350  @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the last iteration is used.
351  @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then 'start_iteration' + 1 is used.
352 
353  Example:
354  @code
355  oscillatory_network = sync(16, type_conn = conn_type.ALL_TO_ALL);
356  output_dynamic = oscillatory_network.simulate_static(100, 10);
357 
358  print("Order parameter at the last step: ", output_dynamic.calculate_order_parameter());
359  print("Order parameter at the first step:", output_dynamic.calculate_order_parameter(0));
360  print("Order parameter evolution between 40 and 50 steps:", output_dynamic.calculate_order_parameter(40, 50));
361  @endcode
362 
363  @return (list) List of levels of global synchronization (order parameter evolution).
364 
365  @see order_estimator
366 
367  """
368 
369  (start_iteration, stop_iteration) = self.__get_start_stop_iterations(start_iteration, stop_iteration);
370 
371  if (self._ccore_sync_dynamic_pointer is not None):
372  return wrapper.sync_dynamic_calculate_order(self._ccore_sync_dynamic_pointer, start_iteration, stop_iteration);
373 
374  sequence_order = [];
375  for index in range(start_iteration, stop_iteration):
376  sequence_order.append(order_estimator.calculate_sync_order(self.output[index]));
377 
378  return sequence_order;
379 
380 
381  def calculate_local_order_parameter(self, oscillatory_network, start_iteration = None, stop_iteration = None):
382  """!
383  @brief Calculates local order parameter.
384  @details Local order parameter or so-called level of local or partial synchronization is calculated by following expression:
385 
386  \f[
387  r_{c}=\left | \sum_{i=0}^{N} \frac{1}{N_{i}} \sum_{j=0}e^{ \theta_{j} - \theta_{i} } \right |;
388  \f]
389 
390  where N - total amount of oscillators in the network and \f$N_{i}\f$ - amount of neighbors of oscillator with index \f$i\f$.
391 
392  @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
393  @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the last iteration is used.
394  @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then 'start_iteration' + 1 is used.
395 
396  @return (list) List of levels of local (partial) synchronization (local order parameter evolution).
397 
398  """
399 
400  (start_iteration, stop_iteration) = self.__get_start_stop_iterations(start_iteration, stop_iteration);
401 
402  if (self._ccore_sync_dynamic_pointer is not None):
403  network_pointer = oscillatory_network._ccore_network_pointer;
404  return wrapper.sync_dynamic_calculate_local_order(self._ccore_sync_dynamic_pointer, network_pointer, start_iteration, stop_iteration);
405 
406  sequence_local_order = [];
407  for index in range(start_iteration, stop_iteration):
408  sequence_local_order.append(order_estimator.calculate_local_sync_order(self.output[index], oscillatory_network));
409 
410  return sequence_local_order;
411 
412 
413  def __get_start_stop_iterations(self, start_iteration, stop_iteration):
414  """!
415  @brief Aplly rules for start_iteration and stop_iteration parameters.
416 
417  @param[in] start_iteration (uint): The first iteration that is used for calculation.
418  @param[in] stop_iteration (uint): The last iteration that is used for calculation.
419 
420  @return (tuple) New the first iteration and the last.
421 
422  """
423  if (start_iteration is None):
424  start_iteration = len(self) - 1;
425 
426  if (stop_iteration is None):
427  stop_iteration = start_iteration + 1;
428 
429  return (start_iteration, stop_iteration);
430 
431 
432 
434  """!
435  @brief Visualizer of output dynamic of sync network (Sync).
436 
437  """
438 
439  @staticmethod
440  def show_output_dynamic(sync_output_dynamic):
441  """!
442  @brief Shows output dynamic (output of each oscillator) during simulation.
443 
444  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
445 
446  @see show_output_dynamics
447 
448  """
449 
450  draw_dynamics(sync_output_dynamic.time, sync_output_dynamic.output, x_title = "t", y_title = "phase", y_lim = [0, 2 * 3.14]);
451 
452 
453  @staticmethod
454  def show_output_dynamics(sync_output_dynamics):
455  """!
456  @brief Shows several output dynamics (output of each oscillator) during simulation.
457  @details Each dynamic is presented on separate plot.
458 
459  @param[in] sync_output_dynamics (list): list of output dynamics 'sync_dynamic' of the Sync network.
460 
461  @see show_output_dynamic
462 
463  """
464 
465  draw_dynamics_set(sync_output_dynamics, "t", "phase", None, [0, 2 * 3.14], False, False);
466 
467 
468  @staticmethod
469  def show_correlation_matrix(sync_output_dynamic, iteration = None):
470  """!
471  @brief Shows correlation matrix between oscillators at the specified iteration.
472 
473  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
474  @param[in] iteration (uint): Number of interation of simulation for which correlation matrix should be allocated.
475  If iternation number is not specified, the last step of simulation is used for the matrix allocation.
476 
477  """
478 
479  _ = plt.figure();
480  correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(iteration);
481 
482  plt.imshow(correlation_matrix, cmap = plt.get_cmap('cool'), interpolation='kaiser', vmin = 0.0, vmax = 1.0);
483  plt.show();
484 
485 
486  @staticmethod
487  def show_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, iteration = None):
488  """!
489  @brief Shows 2D matrix of phase values of oscillators at the specified iteration.
490  @details User should ensure correct matrix sizes in line with following expression grid_width x grid_height that should be equal to
491  amount of oscillators otherwise exception is thrown. If grid_width or grid_height are not specified than phase matrix size
492  will by calculated automatically by square root.
493 
494  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose phase matrix should be shown.
495  @param[in] grid_width (uint): Width of the phase matrix.
496  @param[in] grid_height (uint): Height of the phase matrix.
497  @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
498  If iternation number is not specified, the last step of simulation is used for the matrix allocation.
499 
500  """
501 
502  _ = plt.figure();
503  phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, iteration);
504 
505  plt.imshow(phase_matrix, cmap = plt.get_cmap('jet'), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi);
506  plt.show();
507 
508 
509  @staticmethod
510  def show_order_parameter(sync_output_dynamic, start_iteration = None, stop_iteration = None):
511  """!
512  @brief Shows evolution of order parameter (level of global synchronization in the network).
513 
514  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
515  @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
516  @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
517 
518  """
519 
520  (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
521 
522  order_parameter = sync_output_dynamic.calculate_order_parameter(start_iteration, stop_iteration);
523  axis = plt.subplot(111);
524  plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
525  set_ax_param(axis, "t", "R (order parameter)", None, [0.0, 1.05]);
526 
527  plt.show();
528 
529 
530  @staticmethod
531  def show_local_order_parameter(sync_output_dynamic, oscillatory_network, start_iteration = None, stop_iteration = None):
532  """!
533  @brief Shows evolution of local order parameter (level of local synchronization in the network).
534 
535  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
536  @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
537  @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
538  @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
539 
540  """
541  (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
542 
543  order_parameter = sync_output_dynamic.calculate_local_order_parameter(oscillatory_network, start_iteration, stop_iteration);
544  axis = plt.subplot(111);
545  plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
546  set_ax_param(axis, "t", "R (local order parameter)", None, [0.0, 1.05]);
547 
548  plt.show();
549 
550 
551  @staticmethod
552  def animate_output_dynamic(sync_output_dynamic, animation_velocity = 75, save_movie = None):
553  """!
554  @brief Shows animation of output dynamic (output of each oscillator) during simulation on a circle from [0; 2pi].
555 
556  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
557  @param[in] animation_velocity (uint): Interval between frames in milliseconds.
558  @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
559 
560  """
561 
562  figure = plt.figure();
563 
564  dynamic = sync_output_dynamic.output[0];
565  artist, = plt.polar(dynamic, [1.0] * len(dynamic), 'o', color = 'blue');
566 
567  def init_frame():
568  return [ artist ];
569 
570  def frame_generation(index_dynamic):
571  dynamic = sync_output_dynamic.output[index_dynamic];
572  artist.set_data(dynamic, [1.0] * len(dynamic));
573 
574  return [ artist ];
575 
576  phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = animation_velocity, init_func = init_frame, repeat_delay = 5000);
577 
578  if (save_movie is not None):
579  phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
580  else:
581  plt.show();
582 
583 
584  @staticmethod
585  def animate_correlation_matrix(sync_output_dynamic, animation_velocity = 75, colormap = 'cool', save_movie = None):
586  """!
587  @brief Shows animation of correlation matrix between oscillators during simulation.
588 
589  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
590  @param[in] animation_velocity (uint): Interval between frames in milliseconds.
591  @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
592  @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
593 
594  """
595 
596  figure = plt.figure()
597 
598  correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0)
599  artist = plt.imshow(correlation_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', vmin = 0.0, vmax = 1.0)
600 
601  def init_frame():
602  return [ artist ]
603 
604  def frame_generation(index_dynamic):
605  correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic)
606  artist.set_data(correlation_matrix)
607 
608  return [ artist ]
609 
610  correlation_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000, blit = True)
611 
612  if (save_movie is not None):
613  correlation_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500)
614  else:
615  plt.show()
616 
617 
618  @staticmethod
619  def animate_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, animation_velocity = 75, colormap = 'jet', save_movie = None):
620  """!
621  @brief Shows animation of phase matrix between oscillators during simulation on 2D stage.
622  @details If grid_width or grid_height are not specified than phase matrix size will by calculated automatically by square root.
623 
624  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
625  @param[in] grid_width (uint): Width of the phase matrix.
626  @param[in] grid_height (uint): Height of the phase matrix.
627  @param[in] animation_velocity (uint): Interval between frames in milliseconds.
628  @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
629  @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
630 
631  """
632 
633  figure = plt.figure();
634 
635  def init_frame():
636  return frame_generation(0);
637 
638  def frame_generation(index_dynamic):
639  figure.clf();
640  axis = figure.add_subplot(111);
641 
642  phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, index_dynamic);
643  axis.imshow(phase_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi);
644  artist = figure.gca();
645 
646  return [ artist ];
647 
648  phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000);
649 
650  if (save_movie is not None):
651  phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
652  else:
653  plt.show();
654 
655 
656  @staticmethod
657  def __get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration):
658  """!
659  @brief Apply rule of preparation for start iteration and stop iteration values.
660 
661  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
662  @param[in] start_iteration (uint): The first iteration that is used for calculation.
663  @param[in] stop_iteration (uint): The last iteration that is used for calculation.
664 
665  @return (tuple) New values of start and stop iterations.
666 
667  """
668  if (start_iteration is None):
669  start_iteration = 0;
670 
671  if (stop_iteration is None):
672  stop_iteration = len(sync_output_dynamic);
673 
674  return (start_iteration, stop_iteration);
675 
676 
677  @staticmethod
678  def animate(sync_output_dynamic, title = None, save_movie = None):
679  """!
680  @brief Shows animation of phase coordinates and animation of correlation matrix together for the Sync dynamic output on the same figure.
681 
682  @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
683  @param[in] title (string): Title of the animation that is displayed on a figure if it is specified.
684  @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
685 
686  """
687 
688  dynamic = sync_output_dynamic.output[0];
689  correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0);
690 
691  figure = plt.figure(1);
692  if (title is not None):
693  figure.suptitle(title, fontsize = 26, fontweight = 'bold')
694 
695  ax1 = figure.add_subplot(121, projection='polar');
696  ax2 = figure.add_subplot(122);
697 
698  artist1, = ax1.plot(dynamic, [1.0] * len(dynamic), marker = 'o', color = 'blue', ls = '');
699  artist2 = ax2.imshow(correlation_matrix, cmap = plt.get_cmap('Accent'), interpolation='kaiser');
700 
701  def init_frame():
702  return [ artist1, artist2 ];
703 
704  def frame_generation(index_dynamic):
705  dynamic = sync_output_dynamic.output[index_dynamic];
706  artist1.set_data(dynamic, [1.0] * len(dynamic));
707 
708  correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic);
709  artist2.set_data(correlation_matrix);
710 
711  return [ artist1, artist2 ];
712 
713  dynamic_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = 75, init_func = init_frame, repeat_delay = 5000);
714 
715  if (save_movie is not None):
716  dynamic_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
717  else:
718  plt.show();
719 
720 
721 
723  """!
724  @brief Model of oscillatory network that is based on the Kuramoto model of synchronization.
725 
726  @details CCORE option can be used to use the pyclustering core - C/C++ shared library for processing that significantly increases performance.
727 
728  """
729 
730  def __init__(self, num_osc, weight = 1, frequency = 0, type_conn = conn_type.ALL_TO_ALL, representation = conn_represent.MATRIX, initial_phases = initial_type.RANDOM_GAUSSIAN, ccore = True):
731  """!
732  @brief Constructor of oscillatory network is based on Kuramoto model.
733 
734  @param[in] num_osc (uint): Number of oscillators in the network.
735  @param[in] weight (double): Coupling strength of the links between oscillators.
736  @param[in] frequency (double): Multiplier of internal frequency of the oscillators.
737  @param[in] type_conn (conn_type): Type of connection between oscillators in the network (all-to-all, grid, bidirectional list, etc.).
738  @param[in] representation (conn_represent): Internal representation of connection in the network: matrix or list.
739  @param[in] initial_phases (initial_type): Type of initialization of initial phases of oscillators (random, uniformly distributed, etc.).
740  @param[in] ccore (bool): If True simulation is performed by CCORE library (C++ implementation of pyclustering).
741 
742  """
743 
744  self._ccore_network_pointer = None; # Pointer to CCORE Sync implementation of the network.
745 
746  if ( (ccore is True) and ccore_library.workable() ):
747  self._ccore_network_pointer = wrapper.sync_create_network(num_osc, weight, frequency, type_conn, initial_phases);
748  self._num_osc = num_osc;
749  self._conn_represent = conn_represent.MATRIX;
750 
751  else:
752  super().__init__(num_osc, type_conn, representation);
753 
754  self._weight = weight;
755 
756  self._phases = list();
757  self._freq = list();
758 
759  random.seed();
760  for index in range(0, num_osc, 1):
761  if (initial_phases == initial_type.RANDOM_GAUSSIAN):
762  self._phases.append(random.random() * 2 * pi);
763 
764  elif (initial_phases == initial_type.EQUIPARTITION):
765  self._phases.append( pi / num_osc * index);
766 
767  self._freq.append(random.random() * frequency);
768 
769 
770  def __del__(self):
771  """!
772  @brief Destructor of oscillatory network is based on Kuramoto model.
773 
774  """
775 
776  if (self._ccore_network_pointer is not None):
777  wrapper.sync_destroy_network(self._ccore_network_pointer);
778  self._ccore_network_pointer = None;
779 
780 
781  def sync_order(self):
782  """!
783  @brief Calculates current level of global synchorization (order parameter) in the network.
784  @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when
785  desynchronization is observed in the network. Order parameter is calculated using following equation:
786 
787  \f[
788  r_{c}=\frac{1}{Ne^{i\varphi }}\sum_{j=0}^{N}e^{i\theta_{j}};
789  \f]
790 
791  where \f$\varphi\f$ is a average phase coordinate in the network, \f$N\f$ is an amount of oscillators in the network.
792 
793  Example:
794  @code
795  oscillatory_network = sync(16, type_conn = conn_type.ALL_TO_ALL);
796  output_dynamic = oscillatory_network.simulate_static(100, 10);
797 
798  if (oscillatory_network.sync_order() < 0.9): print("Global synchronization is not reached yet.");
799  else: print("Global synchronization is reached.");
800  @endcode
801 
802  @return (double) Level of global synchronization (order parameter).
803 
804  @see sync_local_order()
805 
806  """
807 
808  if (self._ccore_network_pointer is not None):
809  return wrapper.sync_order(self._ccore_network_pointer);
810 
811  return order_estimator.calculate_sync_order(self._phases);
812 
813 
814  def sync_local_order(self):
815  """!
816  @brief Calculates current level of local (partial) synchronization in the network.
817 
818  @return (double) Level of local (partial) synchronization.
819 
820  @see sync_order()
821 
822  """
823 
824  if (self._ccore_network_pointer is not None):
825  return wrapper.sync_local_order(self._ccore_network_pointer);
826 
827  return order_estimator.calculate_local_sync_order(self._phases, self);
828 
829 
830  def _phase_kuramoto(self, teta, t, argv):
831  """!
832  @brief Returns result of phase calculation for specified oscillator in the network.
833 
834  @param[in] teta (double): Phase of the oscillator that is differentiated.
835  @param[in] t (double): Current time of simulation.
836  @param[in] argv (tuple): Index of the oscillator in the list.
837 
838  @return (double) New phase for specified oscillator (don't assign here).
839 
840  """
841 
842  index = argv;
843  phase = 0;
844  for k in range(0, self._num_osc):
845  if (self.has_connection(index, k) == True):
846  phase += math.sin(self._phases[k] - teta);
847 
848  return ( self._freq[index] + (phase * self._weight / self._num_osc) );
849 
850 
851  def simulate(self, steps, time, solution = solve_type.FAST, collect_dynamic = True):
852  """!
853  @brief Performs static simulation of Sync oscillatory network.
854 
855  @param[in] steps (uint): Number steps of simulations during simulation.
856  @param[in] time (double): Time of simulation.
857  @param[in] solution (solve_type): Type of solution (solving).
858  @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
859 
860  @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
861  otherwise returns only last values (last step of simulation) of dynamic.
862 
863  @see simulate_dynamic()
864  @see simulate_static()
865 
866  """
867 
868  return self.simulate_static(steps, time, solution, collect_dynamic);
869 
870 
871  def simulate_dynamic(self, order = 0.998, solution = solve_type.FAST, collect_dynamic = False, step = 0.1, int_step = 0.01, threshold_changes = 0.0000001):
872  """!
873  @brief Performs dynamic simulation of the network until stop condition is not reached. Stop condition is defined by input argument 'order'.
874 
875  @param[in] order (double): Order of process synchronization, distributed 0..1.
876  @param[in] solution (solve_type): Type of solution.
877  @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
878  @param[in] step (double): Time step of one iteration of simulation.
879  @param[in] int_step (double): Integration step, should be less than step.
880  @param[in] threshold_changes (double): Additional stop condition that helps prevent infinite simulation, defines limit of changes of oscillators between current and previous steps.
881 
882  @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
883  otherwise returns only last values (last step of simulation) of dynamic.
884 
885  @see simulate()
886  @see simulate_static()
887 
888  """
889 
890  if (self._ccore_network_pointer is not None):
891  ccore_instance_dynamic = wrapper.sync_simulate_dynamic(self._ccore_network_pointer, order, solution, collect_dynamic, step, int_step, threshold_changes);
892  return sync_dynamic(None, None, ccore_instance_dynamic);
893 
894  # For statistics and integration
895  time_counter = 0;
896 
897  # Prevent infinite loop. It's possible when required state cannot be reached.
898  previous_order = 0;
899  current_order = self.sync_local_order();
900 
901  # If requested input dynamics
902  dyn_phase = [];
903  dyn_time = [];
904  if (collect_dynamic == True):
905  dyn_phase.append(self._phases);
906  dyn_time.append(0);
907 
908  # Execute until sync state will be reached
909  while (current_order < order):
910  # update states of oscillators
911  self._phases = self._calculate_phases(solution, time_counter, step, int_step);
912 
913  # update time
914  time_counter += step;
915 
916  # if requested input dynamic
917  if (collect_dynamic == True):
918  dyn_phase.append(self._phases);
919  dyn_time.append(time_counter);
920 
921  # update orders
922  previous_order = current_order;
923  current_order = self.sync_local_order();
924 
925  # hang prevention
926  if (abs(current_order - previous_order) < threshold_changes):
927  # print("Warning: sync_network::simulate_dynamic - simulation is aborted due to low level of convergence rate (order = " + str(current_order) + ").");
928  break;
929 
930  if (collect_dynamic != True):
931  dyn_phase.append(self._phases);
932  dyn_time.append(time_counter);
933 
934  output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time, None);
935  return output_sync_dynamic;
936 
937 
938  def simulate_static(self, steps, time, solution = solve_type.FAST, collect_dynamic = False):
939  """!
940  @brief Performs static simulation of oscillatory network.
941 
942  @param[in] steps (uint): Number steps of simulations during simulation.
943  @param[in] time (double): Time of simulation.
944  @param[in] solution (solve_type): Type of solution.
945  @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
946 
947  @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
948  otherwise returns only last values (last step of simulation) of dynamic.
949 
950  @see simulate()
951  @see simulate_dynamic()
952 
953  """
954 
955  if (self._ccore_network_pointer is not None):
956  ccore_instance_dynamic = wrapper.sync_simulate_static(self._ccore_network_pointer, steps, time, solution, collect_dynamic);
957  return sync_dynamic(None, None, ccore_instance_dynamic);
958 
959  dyn_phase = [];
960  dyn_time = [];
961 
962  if (collect_dynamic == True):
963  dyn_phase.append(self._phases);
964  dyn_time.append(0);
965 
966  step = time / steps;
967  int_step = step / 10.0;
968 
969  for t in numpy.arange(step, time + step, step):
970  # update states of oscillators
971  self._phases = self._calculate_phases(solution, t, step, int_step);
972 
973  # update states of oscillators
974  if (collect_dynamic == True):
975  dyn_phase.append(self._phases);
976  dyn_time.append(t);
977 
978  if (collect_dynamic != True):
979  dyn_phase.append(self._phases);
980  dyn_time.append(time);
981 
982  output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time);
983  return output_sync_dynamic;
984 
985 
986  def _calculate_phases(self, solution, t, step, int_step):
987  """!
988  @brief Calculates new phases for oscillators in the network in line with current step.
989 
990  @param[in] solution (solve_type): Type solver of the differential equation.
991  @param[in] t (double): Time of simulation.
992  @param[in] step (double): Step of solution at the end of which states of oscillators should be calculated.
993  @param[in] int_step (double): Step differentiation that is used for solving differential equation.
994 
995  @return (list) New states (phases) for oscillators.
996 
997  """
998 
999  next_phases = [0.0] * self._num_osc; # new oscillator _phases
1000 
1001  for index in range (0, self._num_osc, 1):
1002  if (solution == solve_type.FAST):
1003  result = self._phases[index] + self._phase_kuramoto(self._phases[index], 0, index);
1004  next_phases[index] = self._phase_normalization(result);
1005 
1006  elif ( (solution == solve_type.RK4) or (solution == solve_type.RKF45) ):
1007  result = odeint(self._phase_kuramoto, self._phases[index], numpy.arange(t - step, t, int_step), (index , ));
1008  next_phases[index] = self._phase_normalization(result[len(result) - 1][0]);
1009 
1010  else:
1011  raise NameError("Solver '" + str(solution) + "' is not supported");
1012 
1013  return next_phases;
1014 
1015 
1016  def _phase_normalization(self, teta):
1017  """!
1018  @brief Normalization of phase of oscillator that should be placed between [0; 2 * pi].
1019 
1020  @param[in] teta (double): phase of oscillator.
1021 
1022  @return (double) Normalized phase.
1023 
1024  """
1025 
1026  norm_teta = teta;
1027  while (norm_teta > (2.0 * pi)) or (norm_teta < 0):
1028  if (norm_teta > (2.0 * pi)):
1029  norm_teta -= 2.0 * pi;
1030  else:
1031  norm_teta += 2.0 * pi;
1032 
1033  return norm_teta;
1034 
1035 
1036  def get_neighbors(self, index):
1037  """!
1038  @brief Finds neighbors of the oscillator with specified index.
1039 
1040  @param[in] index (uint): index of oscillator for which neighbors should be found in the network.
1041 
1042  @return (list) Indexes of neighbors of the specified oscillator.
1043 
1044  """
1045 
1046  if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1047  self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1048 
1049  return super().get_neighbors(index);
1050 
1051 
1052  def has_connection(self, i, j):
1053  """!
1054  @brief Returns True if there is connection between i and j oscillators and False - if connection doesn't exist.
1055 
1056  @param[in] i (uint): index of an oscillator in the network.
1057  @param[in] j (uint): index of an oscillator in the network.
1058 
1059  """
1060 
1061  if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1062  self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1063 
1064  return super().has_connection(i, j);
def simulate_dynamic(self, order=0.998, solution=solve_type.FAST, collect_dynamic=False, step=0.1, int_step=0.01, threshold_changes=0.0000001)
Performs dynamic simulation of the network until stop condition is not reached.
Definition: sync.py:871
def animate_correlation_matrix(sync_output_dynamic, animation_velocity=75, colormap='cool', save_movie=None)
Shows animation of correlation matrix between oscillators during simulation.
Definition: sync.py:585
def __del__(self)
Default destructor of Sync dynamic.
Definition: sync.py:158
def __get_start_stop_iterations(self, start_iteration, stop_iteration)
Aplly rules for start_iteration and stop_iteration parameters.
Definition: sync.py:413
def _calculate_phases(self, solution, t, step, int_step)
Calculates new phases for oscillators in the network in line with current step.
Definition: sync.py:986
def calculate_local_sync_order(oscillator_phases, oscillatory_network)
Calculates level of local synchorization (local order parameter) for input phases for the specified n...
Definition: sync.py:84
def time(self)
(list) Returns sampling times when dynamic is measured during simulation.
Definition: sync.py:132
def animate_output_dynamic(sync_output_dynamic, animation_velocity=75, save_movie=None)
Shows animation of output dynamic (output of each oscillator) during simulation on a circle from [0; ...
Definition: sync.py:552
Represents output dynamic of Sync.
Definition: sync.py:113
def show_phase_matrix(sync_output_dynamic, grid_width=None, grid_height=None, iteration=None)
Shows 2D matrix of phase values of oscillators at the specified iteration.
Definition: sync.py:487
def allocate_correlation_matrix(self, iteration=None)
Allocate correlation matrix between oscillators at the specified step of simulation.
Definition: sync.py:302
def simulate_static(self, steps, time, solution=solve_type.FAST, collect_dynamic=False)
Performs static simulation of oscillatory network.
Definition: sync.py:938
Utils that are used by modules of pyclustering.
Definition: __init__.py:1
def sync_local_order(self)
Calculates current level of local (partial) synchronization in the network.
Definition: sync.py:814
def animate_phase_matrix(sync_output_dynamic, grid_width=None, grid_height=None, animation_velocity=75, colormap='jet', save_movie=None)
Shows animation of phase matrix between oscillators during simulation on 2D stage.
Definition: sync.py:619
def allocate_sync_ensembles(self, tolerance=0.01, indexes=None, iteration=None)
Allocate clusters in line with ensembles of synchronous oscillators where each synchronous ensemble c...
Definition: sync.py:194
def show_output_dynamic(sync_output_dynamic)
Shows output dynamic (output of each oscillator) during simulation.
Definition: sync.py:440
def __len__(self)
Returns number of simulation steps that are stored in dynamic.
Definition: sync.py:167
def has_connection(self, i, j)
Returns True if there is connection between i and j oscillators and False - if connection doesn&#39;t exi...
Definition: __init__.py:366
def show_output_dynamics(sync_output_dynamics)
Shows several output dynamics (output of each oscillator) during simulation.
Definition: sync.py:454
def _phase_kuramoto(self, teta, t, argv)
Returns result of phase calculation for specified oscillator in the network.
Definition: sync.py:830
def __getitem__(self, index)
Indexing of the dynamic.
Definition: sync.py:179
def calculate_local_order_parameter(self, oscillatory_network, start_iteration=None, stop_iteration=None)
Calculates local order parameter.
Definition: sync.py:381
def show_correlation_matrix(sync_output_dynamic, iteration=None)
Shows correlation matrix between oscillators at the specified iteration.
Definition: sync.py:469
def calculate_sync_order(oscillator_phases)
Calculates level of global synchronization (order parameter) for input phases.
Definition: sync.py:56
Common network description that consists of information about oscillators and connection between them...
Definition: __init__.py:97
Provides services to calculate order parameter and local order parameter that are used for synchroniz...
Definition: sync.py:49
def __init__(self, phase, time, ccore=None)
Constructor of Sync dynamic.
Definition: sync.py:143
def show_order_parameter(sync_output_dynamic, start_iteration=None, stop_iteration=None)
Shows evolution of order parameter (level of global synchronization in the network).
Definition: sync.py:510
def __del__(self)
Destructor of oscillatory network is based on Kuramoto model.
Definition: sync.py:770
def get_neighbors(self, index)
Finds neighbors of the oscillator with specified index.
Definition: sync.py:1036
def has_connection(self, i, j)
Returns True if there is connection between i and j oscillators and False - if connection doesn&#39;t exi...
Definition: sync.py:1052
def __init__(self, num_osc, weight=1, frequency=0, type_conn=conn_type.ALL_TO_ALL, representation=conn_represent.MATRIX, initial_phases=initial_type.RANDOM_GAUSSIAN, ccore=True)
Constructor of oscillatory network is based on Kuramoto model.
Definition: sync.py:730
def animate(sync_output_dynamic, title=None, save_movie=None)
Shows animation of phase coordinates and animation of correlation matrix together for the Sync dynami...
Definition: sync.py:678
def _phase_normalization(self, teta)
Normalization of phase of oscillator that should be placed between [0; 2 * pi].
Definition: sync.py:1016
def sync_order(self)
Calculates current level of global synchorization (order parameter) in the network.
Definition: sync.py:781
def allocate_phase_matrix(self, grid_width=None, grid_height=None, iteration=None)
Returns 2D matrix of phase values of oscillators at the specified iteration of simulation.
Definition: sync.py:258
def calculate_order_parameter(self, start_iteration=None, stop_iteration=None)
Calculates level of global synchorization (order parameter).
Definition: sync.py:338
Visualizer of output dynamic of sync network (Sync).
Definition: sync.py:433
def show_local_order_parameter(sync_output_dynamic, oscillatory_network, start_iteration=None, stop_iteration=None)
Shows evolution of local order parameter (level of local synchronization in the network).
Definition: sync.py:531
def simulate(self, steps, time, solution=solve_type.FAST, collect_dynamic=True)
Performs static simulation of Sync oscillatory network.
Definition: sync.py:851
def output(self)
(list) Returns output dynamic of the Sync network (phase coordinates of each oscillator in the networ...
Definition: sync.py:120
Neural and oscillatory network module.
Definition: __init__.py:1
Model of oscillatory network that is based on the Kuramoto model of synchronization.
Definition: sync.py:722