Completed
Push — master ( b7ad63...18ed0c )
by Andrei
01:27
created

sync_network.sync_local_order()   A

Complexity

Conditions 2

Size

Total Lines 14

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 2
dl 0
loc 14
rs 9.4285
c 0
b 0
f 0
1
"""!
2
3
@brief Neural Network: Oscillatory Neural Network based on Kuramoto model
4
@details Based on article description:
5
         - A.Arenas, Y.Moreno, C.Zhou. Synchronization in complex networks. 2008.
6
         - X.B.Lu. Adaptive Cluster Synchronization in Coupled Phase Oscillators. 2009.
7
         - X.Lou. Adaptive Synchronizability of Coupled Oscillators With Switching. 2012.
8
         - A.Novikov, E.Benderskaya. Oscillatory Neural Networks Based on the Kuramoto Model. 2014.
9
10
@authors Andrei Novikov ([email protected])
11
@date 2014-2017
12
@copyright GNU Public License
13
14
@cond GNU_PUBLIC_LICENSE
15
    PyClustering is free software: you can redistribute it and/or modify
16
    it under the terms of the GNU General Public License as published by
17
    the Free Software Foundation, either version 3 of the License, or
18
    (at your option) any later version.
19
    
20
    PyClustering is distributed in the hope that it will be useful,
21
    but WITHOUT ANY WARRANTY; without even the implied warranty of
22
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23
    GNU General Public License for more details.
24
    
25
    You should have received a copy of the GNU General Public License
26
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
27
@endcond
28
29
"""
30
31
import matplotlib.pyplot as plt;
0 ignored issues
show
Configuration introduced by
The import matplotlib.pyplot could not be resolved.

This can be caused by one of the following:

1. Missing Dependencies

This error could indicate a configuration issue of Pylint. Make sure that your libraries are available by adding the necessary commands.

# .scrutinizer.yml
before_commands:
    - sudo pip install abc # Python2
    - sudo pip3 install abc # Python3
Tip: We are currently not using virtualenv to run pylint, when installing your modules make sure to use the command for the correct version.

2. Missing __init__.py files

This error could also result from missing __init__.py files in your module folders. Make sure that you place one file in each sub-folder.

Loading history...
32
import matplotlib.animation as animation;
0 ignored issues
show
Configuration introduced by
The import matplotlib.animation could not be resolved.

This can be caused by one of the following:

1. Missing Dependencies

This error could indicate a configuration issue of Pylint. Make sure that your libraries are available by adding the necessary commands.

# .scrutinizer.yml
before_commands:
    - sudo pip install abc # Python2
    - sudo pip3 install abc # Python3
Tip: We are currently not using virtualenv to run pylint, when installing your modules make sure to use the command for the correct version.

2. Missing __init__.py files

This error could also result from missing __init__.py files in your module folders. Make sure that you place one file in each sub-folder.

Loading history...
33
34
import math;
35
import numpy;
0 ignored issues
show
Configuration introduced by
The import numpy could not be resolved.

This can be caused by one of the following:

1. Missing Dependencies

This error could indicate a configuration issue of Pylint. Make sure that your libraries are available by adding the necessary commands.

# .scrutinizer.yml
before_commands:
    - sudo pip install abc # Python2
    - sudo pip3 install abc # Python3
Tip: We are currently not using virtualenv to run pylint, when installing your modules make sure to use the command for the correct version.

2. Missing __init__.py files

This error could also result from missing __init__.py files in your module folders. Make sure that you place one file in each sub-folder.

Loading history...
36
import random;
37
38
import pyclustering.core.sync_wrapper as wrapper;
39
40
from scipy.integrate import odeint;
0 ignored issues
show
Configuration introduced by
The import scipy.integrate could not be resolved.

This can be caused by one of the following:

1. Missing Dependencies

This error could indicate a configuration issue of Pylint. Make sure that your libraries are available by adding the necessary commands.

# .scrutinizer.yml
before_commands:
    - sudo pip install abc # Python2
    - sudo pip3 install abc # Python3
Tip: We are currently not using virtualenv to run pylint, when installing your modules make sure to use the command for the correct version.

2. Missing __init__.py files

This error could also result from missing __init__.py files in your module folders. Make sure that you place one file in each sub-folder.

Loading history...
41
42
from pyclustering.nnet import network, conn_represent, conn_type, initial_type, solve_type;
43
from pyclustering.utils import pi, draw_dynamics, draw_dynamics_set, set_ax_param;
44
45
46
class order_estimator:
47
    """!
48
    @brief Provides services to calculate order parameter and local order parameter that are used
49
            for synchronization level estimation.
50
51
    """
52
    @staticmethod
53
    def calculate_sync_order(oscillator_phases):
54
        """!
55
        @brief Calculates level of global synchronization (order parameter) for input phases.
56
        @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when 
57
                  desynchronization is observed in the network.
58
        
59
        @param[in] oscillator_phases (list): List of oscillator phases that are used for level of global synchronization.
60
        
61
        @return (double) Level of global synchronization (order parameter).
62
        
63
        @see calculate_order_parameter()
64
        
65
        """
66
        
67
        exp_amount = 0.0;
68
        average_phase = 0.0;
69
70
        for phase in oscillator_phases:
71
            exp_amount += math.expm1( abs(1j * phase) );
72
            average_phase += phase;
73
        
74
        exp_amount /= len(oscillator_phases);
75
        average_phase = math.expm1( abs(1j * (average_phase / len(oscillator_phases))) );
76
        
77
        return abs(average_phase) / abs(exp_amount);
78
79
80
    @staticmethod
81
    def calculate_local_sync_order(oscillator_phases, oscillatory_network):
82
        """!
83
        @brief Calculates level of local synchorization (local order parameter) for input phases for the specified network.
84
        @details This parameter is tend 1.0 when the oscillatory network close to local synchronization and it tend to 0.0 when 
85
                  desynchronization is observed in the network.
86
        
87
        @param[in] oscillator_phases (list): List of oscillator phases that are used for level of local (partial) synchronization.
88
        @param[in] oscillatory_network (sync): Instance of oscillatory network whose connections are required for calculation.
89
        
90
        @return (double) Level of local synchronization (local order parameter).
91
        
92
        """
93
        
94
        exp_amount = 0.0;
95
        num_neigh = 0.0;
96
        
97
        for i in range(0, len(oscillatory_network), 1):
98
            for j in range(0, len(oscillatory_network), 1):
99
                if (oscillatory_network.has_connection(i, j) == True):
100
                    exp_amount += math.exp(-abs(oscillator_phases[j] - oscillator_phases[i]));
101
                    num_neigh += 1.0;
102
        
103
        if (num_neigh == 0):
104
            num_neigh = 1.0;
105
        
106
        return exp_amount / num_neigh;
107
108
109
110
class sync_dynamic:
111
    """!
112
    @brief Represents output dynamic of Sync.
113
    
114
    """
115
    
116
    @property
117
    def output(self):
118
        """!
119
        @brief (list) Returns output dynamic of the Sync network (phase coordinates of each oscillator in the network) during simulation.
120
        
121
        """
122
        if ( (self._ccore_sync_dynamic_pointer is not None) and ( (self._dynamic is None) or (len(self._dynamic) == 0) ) ):
123
            self._dynamic = wrapper.sync_dynamic_get_output(self._ccore_sync_dynamic_pointer);
124
        
125
        return self._dynamic;
126
    
127
    
128
    @property
129
    def time(self):
130
        """!
131
        @brief (list) Returns sampling times when dynamic is measured during simulation.
132
        
133
        """
134
        if ( (self._ccore_sync_dynamic_pointer is not None) and ( (self._time is None) or (len(self._time) == 0) ) ):
135
            self._time = wrapper.sync_dynamic_get_time(self._ccore_sync_dynamic_pointer);
136
        
137
        return self._time;
138
    
139
    
140
    def __init__(self, phase, time, ccore = None):
141
        """!
142
        @brief Constructor of Sync dynamic.
143
        
144
        @param[in] phase (list): Dynamic of oscillators on each step of simulation. If ccore pointer is specified than it can be ignored.
145
        @param[in] time (list): Simulation time.
146
        @param[in] ccore (ctypes.pointer): Pointer to CCORE sync_dynamic instance in memory.
147
        
148
        """
149
        
150
        self._dynamic = phase;
151
        self._time = time;
152
        self._ccore_sync_dynamic_pointer = ccore;
153
    
154
    
155
    def __del__(self):
156
        """!
157
        @brief Default destructor of Sync dynamic.
158
        
159
        """
160
        if (self._ccore_sync_dynamic_pointer is not None):
161
            wrapper.sync_dynamic_destroy(self._ccore_sync_dynamic_pointer);
162
    
163
    
164
    def __len__(self):
165
        """!
166
        @brief (uint) Returns number of simulation steps that are stored in dynamic.
167
        
168
        """
169
        if (self._ccore_sync_dynamic_pointer is not None):
170
            return wrapper.sync_dynamic_get_size(self._ccore_sync_dynamic_pointer);
171
        
172
        return len(self._dynamic);
173
    
174
    
175
    def __getitem__(self, index):
176
        """!
177
        @brief Indexing of the dynamic.
178
        
179
        """
180
        if (index is 0):
181
            return self.time;
182
        
183
        elif (index is 1):
184
            return self.output;
185
        
186
        else:
187
            raise NameError('Out of range ' + index + ': only indexes 0 and 1 are supported.');
188
189
190
    def allocate_sync_ensembles(self, tolerance = 0.01, indexes = None, iteration = None):
191
        """!
192
        @brief Allocate clusters in line with ensembles of synchronous oscillators where each synchronous ensemble corresponds to only one cluster.
193
               
194
        @param[in] tolerance (double): Maximum error for allocation of synchronous ensemble oscillators.
195
        @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]).
196
        @param[in] iteration (uint): Iteration of simulation that should be used for allocation.
197
        
198
        @return (list) Grours (lists) of indexes of synchronous oscillators.
199
                For example [ [index_osc1, index_osc3], [index_osc2], [index_osc4, index_osc5] ].
200
        
201
        """
202
        
203
        if (self._ccore_sync_dynamic_pointer is not None):
204
            ensembles = wrapper.sync_dynamic_allocate_sync_ensembles(self._ccore_sync_dynamic_pointer, tolerance, iteration);
205
            
206
            if (indexes is not None):
207
                for ensemble in ensembles:
208
                    for index in range(len(ensemble)):
209
                        ensemble[index] = indexes[ ensemble[index] ];
210
                
211
            return ensembles;
212
        
213
        if ( (self._dynamic is None) or (len(self._dynamic) == 0) ):
214
            return [];
215
        
216
        number_oscillators = len(self._dynamic[0]);
217
        last_state = None;
218
        
219
        if (iteration is None):
220
            last_state = self._dynamic[len(self._dynamic) - 1];
221
        else:
222
            last_state = self._dynamic[iteration];
223
        
224
        clusters = [];
225
        if (number_oscillators > 0):
226
            clusters.append([0]);
227
        
228
        for i in range(1, number_oscillators, 1):
229
            cluster_allocated = False;
230
            for cluster in clusters:
231
                for neuron_index in cluster:
232
                    last_state_shifted = abs(last_state[i] - 2 * pi);
233
                    
234
                    if ( ( (last_state[i] < (last_state[neuron_index] + tolerance)) and (last_state[i] > (last_state[neuron_index] - tolerance)) ) or
235
                         ( (last_state_shifted < (last_state[neuron_index] + tolerance)) and (last_state_shifted > (last_state[neuron_index] - tolerance)) ) ):
236
                        cluster_allocated = True;
237
                        
238
                        real_index = i;
239
                        if (indexes is not None):
240
                            real_index = indexes[i];
241
                        
242
                        cluster.append(real_index);
243
                        break;
244
                
245
                if (cluster_allocated == True):
246
                    break;
247
            
248
            if (cluster_allocated == False):
249
                clusters.append([i]);
250
        
251
        return clusters;
252
    
253
    
254
    def allocate_phase_matrix(self, grid_width = None, grid_height = None, iteration = None):
255
        """!
256
        @brief Returns 2D matrix of phase values of oscillators at the specified iteration of simulation.
257
        @details User should ensure correct matrix sizes in line with following expression grid_width x grid_height that should be equal to 
258
                  amount of oscillators otherwise exception is thrown. If grid_width or grid_height are not specified than phase matrix size 
259
                  will by calculated automatically by square root.
260
        
261
        @param[in] grid_width (uint): Width of the allocated matrix.
262
        @param[in] grid_height (uint): Height of the allocated matrix.
263
        @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
264
                    If iternation number is not specified, the last step of simulation is used for the matrix allocation.
265
        
266
        @return (list) Phase value matrix of oscillators with size [number_oscillators x number_oscillators].
267
        
268
        """
269
        
270
        output_dynamic = self.output;
271
        
272
        if ( (output_dynamic is None) or (len(output_dynamic) == 0) ):
273
            return [];
274
        
275
        current_dynamic = output_dynamic[len(output_dynamic) - 1];
276
        if (iteration is not None):
277
            current_dynamic = output_dynamic[iteration];
278
        
279
        width_matrix = grid_width;
280
        height_matrix = grid_height;
281
        number_oscillators = len(current_dynamic);
282
        if ( (width_matrix is None) or (height_matrix is None) ):
283
            width_matrix = int(math.ceil(math.sqrt(number_oscillators)));
284
            height_matrix = width_matrix;
285
286
        if (number_oscillators != width_matrix * height_matrix):
287
            raise NameError("Impossible to allocate phase matrix with specified sizes, amout of neurons should be equal to grid_width * grid_height.");
288
        
289
        
290
        phase_matrix = [ [ 0.0 for i in range(width_matrix) ] for j in range(height_matrix) ];
291
        for i in range(height_matrix):
292
            for j in range(width_matrix):
293
                phase_matrix[i][j] = current_dynamic[j + i * width_matrix];
294
        
295
        return phase_matrix;
296
    
297
    
298
    def allocate_correlation_matrix(self, iteration = None):
299
        """!
300
        @brief Allocate correlation matrix between oscillators at the specified step of simulation.
301
               
302
        @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
303
                    If iternation number is not specified, the last step of simulation is used for the matrix allocation.
304
        
305
        @return (list) Correlation matrix between oscillators with size [number_oscillators x number_oscillators].
306
        
307
        """
308
        
309
        if (self._ccore_sync_dynamic_pointer is not None):
310
            return wrapper.sync_dynamic_allocate_correlation_matrix(self._ccore_sync_dynamic_pointer, iteration);
311
        
312
        if ( (self._dynamic is None) or (len(self._dynamic) == 0) ):
313
            return [];
314
        
315
        dynamic = self._dynamic;
316
        current_dynamic = dynamic[len(dynamic) - 1];
317
        
318
        if (iteration is not None):
319
            current_dynamic = dynamic[iteration];
320
        
321
        number_oscillators = len(dynamic[0]);
322
        affinity_matrix = [ [ 0.0 for i in range(number_oscillators) ] for j in range(number_oscillators) ];
323
        
324
        for i in range(number_oscillators):
325
            for j in range(number_oscillators):
326
                phase1 = current_dynamic[i];
327
                phase2 = current_dynamic[j];
328
                
329
                affinity_matrix[i][j] = abs(math.sin(phase1 - phase2));
330
                
331
        return affinity_matrix;
332
333
334
    def calculate_order_parameter(self, start_iteration = None, stop_iteration = None):
335
        """!
0 ignored issues
show
Bug introduced by
A suspicious escape sequence \s was found. Did you maybe forget to add an r prefix?

Escape sequences in Python are generally interpreted according to rules similar to standard C. Only if strings are prefixed with r or R are they interpreted as regular expressions.

The escape sequence that was used indicates that you might have intended to write a regular expression.

Learn more about the available escape sequences. in the Python documentation.

Loading history...
336
        @brief Calculates level of global synchorization (order parameter).
337
        @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when 
338
                  desynchronization is observed in the network. Order parameter is calculated using following equation:
339
                  
340
                  \f[
341
                  r_{c}=\frac{1}{Ne^{i\varphi }}\sum_{j=0}^{N}e^{i\theta_{j}};
342
                  \f]
343
                  
344
                  where \f$\varphi\f$ is a average phase coordinate in the network, \f$N\f$ is an amount of oscillators in the network.
345
        
346
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the last iteration is used.
347
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then 'start_iteration' + 1 is used.
348
        
349
        Example:
350
        @code
351
            oscillatory_network = sync(16, type_conn = conn_type.ALL_TO_ALL);
352
            output_dynamic = oscillatory_network.simulate_static(100, 10);
353
            
354
            print("Order parameter at the last step: ", output_dynamic.calculate_order_parameter());
355
            print("Order parameter at the first step:", output_dynamic.calculate_order_parameter(0));
356
            print("Order parameter evolution between 40 and 50 steps:", output_dynamic.calculate_order_parameter(40, 50));
357
        @endcode
358
        
359
        @return (list) List of levels of global synchronization (order parameter evolution).
360
        
361
        @see order_estimator
362
        
363
        """
364
        
365
        (start_iteration, stop_iteration) = self.__get_start_stop_iterations(start_iteration, stop_iteration);
366
        
367
        if (self._ccore_sync_dynamic_pointer is not None):
368
            return wrapper.sync_dynamic_calculate_order(self._ccore_sync_dynamic_pointer, start_iteration, stop_iteration);
369
        
370
        sequence_order = [];
371
        for index in range(start_iteration, stop_iteration):
372
            sequence_order.append(order_estimator.calculate_sync_order(self.output[index]));
373
        
374
        return sequence_order;
375
376
377
    def calculate_local_order_parameter(self, oscillatory_network, start_iteration = None, stop_iteration = None):
378
        """!
0 ignored issues
show
Bug introduced by
A suspicious escape sequence \l was found. Did you maybe forget to add an r prefix?

Escape sequences in Python are generally interpreted according to rules similar to standard C. Only if strings are prefixed with r or R are they interpreted as regular expressions.

The escape sequence that was used indicates that you might have intended to write a regular expression.

Learn more about the available escape sequences. in the Python documentation.

Loading history...
Bug introduced by
A suspicious escape sequence \s was found. Did you maybe forget to add an r prefix?

Escape sequences in Python are generally interpreted according to rules similar to standard C. Only if strings are prefixed with r or R are they interpreted as regular expressions.

The escape sequence that was used indicates that you might have intended to write a regular expression.

Learn more about the available escape sequences. in the Python documentation.

Loading history...
379
        @brief Calculates local order parameter.
380 View Code Duplication
        @details Local order parameter or so-called level of local or partial synchronization is calculated by following expression:
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
381
        
382
        \f[
383
        r_{c}=\left | \sum_{i=0}^{N} \frac{1}{N_{i}} \sum_{j=0}e^{ \theta_{j} - \theta_{i} } \right |;
384
        \f]
385
        
386
        where N - total amount of oscillators in the network and \f$N_{i}\f$ - amount of neighbors of oscillator with index \f$i\f$.
387
        
388
        @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
389
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the last iteration is used.
390
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then 'start_iteration' + 1 is used.
391
        
392
        @return (list) List of levels of local (partial) synchronization (local order parameter evolution).
393
        
394
        """
395
396
        (start_iteration, stop_iteration) = self.__get_start_stop_iterations(start_iteration, stop_iteration);
397
        
398
        if (self._ccore_sync_dynamic_pointer is not None):
399
            network_pointer = oscillatory_network._ccore_network_pointer;
400
            return wrapper.sync_dynamic_calculate_local_order(self._ccore_sync_dynamic_pointer, network_pointer, start_iteration, stop_iteration);
401
        
402
        sequence_local_order = [];
403
        for index in range(start_iteration, stop_iteration):
404
            sequence_local_order.append(order_estimator.calculate_local_sync_order(self.output[index], oscillatory_network));
405
        
406
        return sequence_local_order;
407
408
409
    def __get_start_stop_iterations(self, start_iteration, stop_iteration):
410
        """!
411
        @brief Aplly rules for start_iteration and stop_iteration parameters.
412
413
        @param[in] start_iteration (uint): The first iteration that is used for calculation.
414 View Code Duplication
        @param[in] stop_iteration (uint): The last iteration that is used for calculation.
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
415
        
416
        @return (tuple) New the first iteration and the last.
417
        
418
        """
419
        if (start_iteration is None):
420
            start_iteration = len(self) - 1;
421
        
422
        if (stop_iteration is None):
423
            stop_iteration = start_iteration + 1;
424
        
425
        return (start_iteration, stop_iteration);
426
427
428
429
class sync_visualizer:
430
    """!
431
    @brief Visualizer of output dynamic of sync network (Sync).
432
    
433
    """
434
435
    @staticmethod
436
    def show_output_dynamic(sync_output_dynamic):
437
        """!
438
        @brief Shows output dynamic (output of each oscillator) during simulation.
439
        
440
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
441
        
442
        @see show_output_dynamics
443
        
444
        """
445
        
446
        draw_dynamics(sync_output_dynamic.time, sync_output_dynamic.output, x_title = "t", y_title = "phase", y_lim = [0, 2 * 3.14]);
447
    
448
    
449
    @staticmethod
450
    def show_output_dynamics(sync_output_dynamics):
451
        """!
452
        @brief Shows several output dynamics (output of each oscillator) during simulation.
453
        @details Each dynamic is presented on separate plot.
454
        
455
        @param[in] sync_output_dynamics (list): list of output dynamics 'sync_dynamic' of the Sync network.
456
        
457
        @see show_output_dynamic
458
        
459
        """
460
        
461
        draw_dynamics_set(sync_output_dynamics, "t", "phase", None, [0, 2 * 3.14], False, False);
462
    
463
    
464
    @staticmethod
465
    def show_correlation_matrix(sync_output_dynamic, iteration = None):
466
        """!
467
        @brief Shows correlation matrix between oscillators at the specified iteration.
468
        
469
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
470
        @param[in] iteration (uint): Number of interation of simulation for which correlation matrix should be allocated.
471
                                      If iternation number is not specified, the last step of simulation is used for the matrix allocation.
472
        
473
        """
474
        
475
        _ = plt.figure();
476
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(iteration);
477
        
478
        plt.imshow(correlation_matrix, cmap = plt.get_cmap('cool'), interpolation='kaiser', vmin = 0.0, vmax = 1.0); 
479
        plt.show();
480
481
482
    @staticmethod
483
    def show_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, iteration = None):
484
        """!
485
        @brief Shows 2D matrix of phase values of oscillators at the specified iteration.
486
        @details User should ensure correct matrix sizes in line with following expression grid_width x grid_height that should be equal to 
487
                  amount of oscillators otherwise exception is thrown. If grid_width or grid_height are not specified than phase matrix size 
488
                  will by calculated automatically by square root.
489
        
490
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose phase matrix should be shown.
491
        @param[in] grid_width (uint): Width of the phase matrix.
492
        @param[in] grid_height (uint): Height of the phase matrix.
493
        @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
494
                    If iternation number is not specified, the last step of simulation is used for the matrix allocation.
495
        
496
        """
497
        
498
        _ = plt.figure();
499
        phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, iteration);
500
        
501
        plt.imshow(phase_matrix, cmap = plt.get_cmap('jet'), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi); 
502
        plt.show();
503
504
505
    @staticmethod
506
    def show_order_parameter(sync_output_dynamic, start_iteration = None, stop_iteration = None):
507
        """!
508
        @brief Shows evolution of order parameter (level of global synchronization in the network).
509
        
510
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
511
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
512
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
513
        
514
        """
515
        
516
        (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
517
        
518
        order_parameter = sync_output_dynamic.calculate_order_parameter(start_iteration, stop_iteration);
519
        axis = plt.subplot(111);
520
        plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
521
        set_ax_param(axis, "t", "R (order parameter)", None, [0.0, 1.05]);
522
        
523
        plt.show();
524
525
526
    @staticmethod
527
    def show_local_order_parameter(sync_output_dynamic, oscillatory_network, start_iteration = None, stop_iteration = None):
528
        """!
529
        @brief Shows evolution of local order parameter (level of local synchronization in the network).
530
        
531
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
532
        @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
533
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
534
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
535
        
536
        """
537
        (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
538
        
539
        order_parameter = sync_output_dynamic.calculate_local_order_parameter(oscillatory_network, start_iteration, stop_iteration);
540
        axis = plt.subplot(111);
541
        plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
542
        set_ax_param(axis, "t", "R (local order parameter)", None, [0.0, 1.05]);
543
        
544
        plt.show();
545
546
547
    @staticmethod
548
    def animate_output_dynamic(sync_output_dynamic, animation_velocity = 75, save_movie = None):
549
        """!
550
        @brief Shows animation of output dynamic (output of each oscillator) during simulation on a circle from [0; 2pi].
551
        
552
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
553
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
554
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
555
        
556
        """
557
        
558
        figure = plt.figure();
559
        
560
        dynamic = sync_output_dynamic.output[0];
561
        artist, = plt.polar(dynamic, [1.0] * len(dynamic), 'o', color = 'blue');
562
        
563
        def init_frame():
564
            return [ artist ];
565
        
566
        def frame_generation(index_dynamic):
567
            dynamic = sync_output_dynamic.output[index_dynamic];
568
            artist.set_data(dynamic, [1.0] * len(dynamic));
569
            
570
            return [ artist ];
571
        
572
        phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = animation_velocity, init_func = init_frame, repeat_delay = 5000);
573
574
        if (save_movie is not None):
575
            phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
576
        else:
577
            plt.show();
578
579
580
    @staticmethod
581
    def animate_correlation_matrix(sync_output_dynamic, animation_velocity = 75, colormap = 'cool', save_movie = None):
582
        """!
583
        @brief Shows animation of correlation matrix between oscillators during simulation.
584
        
585
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
586
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
587
        @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
588
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
589
        
590
        """
591
        
592
        figure = plt.figure();
593
        
594
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0);
595
        artist = plt.imshow(correlation_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', hold = True, vmin = 0.0, vmax = 1.0);
596
        
597
        def init_frame(): 
598
            return [ artist ];
599
        
600
        def frame_generation(index_dynamic):
601
            correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic);
602
            artist.set_data(correlation_matrix);
603
            
604
            return [ artist ];
605
606
        correlation_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000, blit = True);
607
        
608
        if (save_movie is not None):
609
            correlation_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
610
        else:
611
            plt.show();
612
613
614
    @staticmethod
615
    def animate_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, animation_velocity = 75, colormap = 'jet', save_movie = None):
616
        """!
617
        @brief Shows animation of phase matrix between oscillators during simulation on 2D stage.
618
        @details If grid_width or grid_height are not specified than phase matrix size will by calculated automatically by square root.
619
        
620
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
621
        @param[in] grid_width (uint): Width of the phase matrix.
622
        @param[in] grid_height (uint): Height of the phase matrix.
623
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
624
        @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
625
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
626
        
627
        """
628
        
629
        figure = plt.figure();
630
        
631
        def init_frame(): 
632
            return frame_generation(0);
633
        
634
        def frame_generation(index_dynamic):
635
            figure.clf();
636
            axis = figure.add_subplot(111);
637
            
638
            phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, index_dynamic);
639
            axis.imshow(phase_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi);
640
            artist = figure.gca();
641
            
642
            return [ artist ];
643
644
        phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000);
645
        
646
        if (save_movie is not None):
647
            phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
648
        else:
649
            plt.show();
650
651
652
    @staticmethod
653
    def __get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration):
654
        """!
655
        @brief Apply rule of preparation for start iteration and stop iteration values.
656
        
657
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
658
        @param[in] start_iteration (uint): The first iteration that is used for calculation.
659
        @param[in] stop_iteration (uint): The last iteration that is used for calculation.
660
        
661
        @return (tuple) New values of start and stop iterations.
662
        
663
        """
664
        if (start_iteration is None):
665
            start_iteration = 0;
666
        
667
        if (stop_iteration is None):
668
            stop_iteration = len(sync_output_dynamic);
669
        
670
        return (start_iteration, stop_iteration);
671
672
673
    @staticmethod
674
    def animate(sync_output_dynamic, title = None, save_movie = None):
675
        """!
676
        @brief Shows animation of phase coordinates and animation of correlation matrix together for the Sync dynamic output on the same figure.
677
        
678
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
679
        @param[in] title (string): Title of the animation that is displayed on a figure if it is specified.
680
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
681
        
682
        """
683
        
684
        dynamic = sync_output_dynamic.output[0];
685
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0);
686
        
687
        figure = plt.figure(1);
688
        if (title is not None):
689
            figure.suptitle(title, fontsize = 26, fontweight = 'bold')
690
        
691
        ax1 = figure.add_subplot(121, projection='polar');
692
        ax2 = figure.add_subplot(122);
693
        
694
        artist1, = ax1.plot(dynamic, [1.0] * len(dynamic), marker = 'o', color = 'blue', ls = '');
695
        artist2 = ax2.imshow(correlation_matrix, cmap = plt.get_cmap('Accent'), interpolation='kaiser');
696
        
697
        def init_frame():
698
            return [ artist1, artist2 ];
699
700
        def frame_generation(index_dynamic):
701
            dynamic = sync_output_dynamic.output[index_dynamic];
702
            artist1.set_data(dynamic, [1.0] * len(dynamic));
703
            
704
            correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic);
705
            artist2.set_data(correlation_matrix);
706
            
707
            return [ artist1, artist2 ];
708
        
709
        dynamic_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = 75, init_func = init_frame, repeat_delay = 5000);
710
        
711
        if (save_movie is not None):
712
            dynamic_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
713
        else:
714
            plt.show();
715
716
717
718
class sync_network(network):
719
    """!
720
    @brief Model of oscillatory network that is based on the Kuramoto model of synchronization.
721
    
722
    """
723
724
    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 = False):
725
        """!
726
        @brief Constructor of oscillatory network is based on Kuramoto model.
727
        
728
        @param[in] num_osc (uint): Number of oscillators in the network.
729
        @param[in] weight (double): Coupling strength of the links between oscillators.
730
        @param[in] frequency (double): Multiplier of internal frequency of the oscillators.
731
        @param[in] type_conn (conn_type): Type of connection between oscillators in the network (all-to-all, grid, bidirectional list, etc.).
732
        @param[in] representation (conn_represent): Internal representation of connection in the network: matrix or list.
733
        @param[in] initial_phases (initial_type): Type of initialization of initial phases of oscillators (random, uniformly distributed, etc.).
734
        @param[in] ccore (bool): If True simulation is performed by CCORE library (C++ implementation of pyclustering).
735
        
736
        """
737
        
738
        self._ccore_network_pointer = None;      # Pointer to CCORE Sync implementation of the network.
739
        
740
        if (ccore is True):
741
            self._ccore_network_pointer = wrapper.sync_create_network(num_osc, weight, frequency, type_conn, initial_phases);
742
            self._num_osc = num_osc;
743
            self._conn_represent = conn_represent.MATRIX;
744
        
745
        else:   
746
            super().__init__(num_osc, type_conn, representation);
747
            
748
            self._weight = weight;
749
            
750
            self._phases = list();
751
            self._freq = list();
752
            
753
            random.seed();
754
            for index in range(0, num_osc, 1):
755
                if (initial_phases == initial_type.RANDOM_GAUSSIAN):
756
                    self._phases.append(random.random() * 2 * pi);
757
                
758
                elif (initial_phases == initial_type.EQUIPARTITION):
759
                    self._phases.append( pi / num_osc * index);
760
                
761
                self._freq.append(random.random() * frequency);
762
763
764
    def __del__(self):
765
        """!
766
        @brief Destructor of oscillatory network is based on Kuramoto model.
767
        
768
        """
769
        
770
        if (self._ccore_network_pointer is not None):
771
            wrapper.sync_destroy_network(self._ccore_network_pointer);
772
            self._ccore_network_pointer = None;
773
774
775
    def sync_order(self):
776
        """!
0 ignored issues
show
Bug introduced by
A suspicious escape sequence \s was found. Did you maybe forget to add an r prefix?

Escape sequences in Python are generally interpreted according to rules similar to standard C. Only if strings are prefixed with r or R are they interpreted as regular expressions.

The escape sequence that was used indicates that you might have intended to write a regular expression.

Learn more about the available escape sequences. in the Python documentation.

Loading history...
777
        @brief Calculates current level of global synchorization (order parameter) in the network.
778
        @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when 
779
                  desynchronization is observed in the network. Order parameter is calculated using following equation:
780
                  
781
                  \f[
782
                  r_{c}=\frac{1}{Ne^{i\varphi }}\sum_{j=0}^{N}e^{i\theta_{j}};
783
                  \f]
784
                  
785
                  where \f$\varphi\f$ is a average phase coordinate in the network, \f$N\f$ is an amount of oscillators in the network.
786
        
787
        Example:
788
        @code
789
            oscillatory_network = sync(16, type_conn = conn_type.ALL_TO_ALL);
790
            output_dynamic = oscillatory_network.simulate_static(100, 10);
791
            
792
            if (oscillatory_network.sync_order() < 0.9): print("Global synchronization is not reached yet.");
793
            else: print("Global synchronization is reached.");
794
        @endcode
795
        
796
        @return (double) Level of global synchronization (order parameter).
797
        
798
        @see sync_local_order()
799
        
800
        """
801
        
802
        if (self._ccore_network_pointer is not None):
803
            return wrapper.sync_order(self._ccore_network_pointer);
804
        
805
        return order_estimator.calculate_sync_order(self._phases);
806
807
808
    def sync_local_order(self):
809
        """!
810
        @brief Calculates current level of local (partial) synchronization in the network.
811
        
812
        @return (double) Level of local (partial) synchronization.
813
        
814
        @see sync_order()
815
        
816
        """
817
        
818
        if (self._ccore_network_pointer is not None):
819
            return wrapper.sync_local_order(self._ccore_network_pointer);
820
        
821
        return order_estimator.calculate_local_sync_order(self._phases, self);
822
823
824
    def _phase_kuramoto(self, teta, t, argv):
0 ignored issues
show
Unused Code introduced by
The argument t seems to be unused.
Loading history...
825
        """!
826
        @brief Returns result of phase calculation for specified oscillator in the network.
827
        
828
        @param[in] teta (double): Phase of the oscillator that is differentiated.
829
        @param[in] t (double): Current time of simulation.
830
        @param[in] argv (tuple): Index of the oscillator in the list.
831
        
832
        @return (double) New phase for specified oscillator (don't assign here).
833
        
834
        """
835
        
836
        index = argv;
837
        phase = 0;
838
        for k in range(0, self._num_osc):
839
            if (self.has_connection(index, k) == True):
840
                phase += math.sin(self._phases[k] - teta);
841
            
842
        return ( self._freq[index] + (phase * self._weight / self._num_osc) );
843
844
845
    def simulate(self, steps, time, solution = solve_type.FAST, collect_dynamic = True):
846
        """!
847
        @brief Performs static simulation of Sync oscillatory network.
848
        
849
        @param[in] steps (uint): Number steps of simulations during simulation.
850
        @param[in] time (double): Time of simulation.
851
        @param[in] solution (solve_type): Type of solution (solving).
852
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
853
        
854
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
855
                otherwise returns only last values (last step of simulation) of dynamic.
856
        
857
        @see simulate_dynamic()
858
        @see simulate_static()
859
        
860
        """
861
        
862
        return self.simulate_static(steps, time, solution, collect_dynamic);
863
864
865
    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):
866
        """!
867
        @brief Performs dynamic simulation of the network until stop condition is not reached. Stop condition is defined by input argument 'order'.
868
        
869
        @param[in] order (double): Order of process synchronization, distributed 0..1.
870
        @param[in] solution (solve_type): Type of solution.
871
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
872
        @param[in] step (double): Time step of one iteration of simulation.
873
        @param[in] int_step (double): Integration step, should be less than step.
874
        @param[in] threshold_changes (double): Additional stop condition that helps prevent infinite simulation, defines limit of changes of oscillators between current and previous steps.
875
        
876
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
877
                otherwise returns only last values (last step of simulation) of dynamic.
878
        
879
        @see simulate()
880
        @see simulate_static()
881
        
882
        """
883
        
884
        if (self._ccore_network_pointer is not None):
885
            ccore_instance_dynamic = wrapper.sync_simulate_dynamic(self._ccore_network_pointer, order, solution, collect_dynamic, step, int_step, threshold_changes);
886
            return sync_dynamic(None, None, ccore_instance_dynamic);
887
        
888
        # For statistics and integration
889
        time_counter = 0;
890
        
891
        # Prevent infinite loop. It's possible when required state cannot be reached.
892
        previous_order = 0;
893
        current_order = self.sync_local_order();
894
        
895
        # If requested input dynamics
896
        dyn_phase = [];
897
        dyn_time = [];
898
        if (collect_dynamic == True):
899
            dyn_phase.append(self._phases);
900
            dyn_time.append(0);
901
        
902
        # Execute until sync state will be reached
903
        while (current_order < order):
904
            # update states of oscillators
905
            self._phases = self._calculate_phases(solution, time_counter, step, int_step);
906
            
907
            # update time
908
            time_counter += step;
909
            
910
            # if requested input dynamic
911
            if (collect_dynamic == True):
912
                dyn_phase.append(self._phases);
913
                dyn_time.append(time_counter);
914
                
915
            # update orders
916
            previous_order = current_order;
917
            current_order = self.sync_local_order();
918
            
919
            # hang prevention
920
            if (abs(current_order - previous_order) < threshold_changes):
921
                # print("Warning: sync_network::simulate_dynamic - simulation is aborted due to low level of convergence rate (order = " + str(current_order) + ").");
922
                break;
923
            
924
        if (collect_dynamic != True):
925
            dyn_phase.append(self._phases);
926
            dyn_time.append(time_counter);
927
928
        output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time, None);
929
        return output_sync_dynamic;
930
931
932
    def simulate_static(self, steps, time, solution = solve_type.FAST, collect_dynamic = False):
933
        """!
934
        @brief Performs static simulation of oscillatory network.
935
        
936
        @param[in] steps (uint): Number steps of simulations during simulation.
937
        @param[in] time (double): Time of simulation.
938
        @param[in] solution (solve_type): Type of solution.
939
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
940
        
941
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
942
                otherwise returns only last values (last step of simulation) of dynamic.
943
        
944
        @see simulate()
945
        @see simulate_dynamic()
946
        
947
        """
948
        
949
        if (self._ccore_network_pointer is not None):
950
            ccore_instance_dynamic = wrapper.sync_simulate_static(self._ccore_network_pointer, steps, time, solution, collect_dynamic);
951
            return sync_dynamic(None, None, ccore_instance_dynamic);
952
        
953
        dyn_phase = [];
954
        dyn_time = [];
955
        
956
        if (collect_dynamic == True):
957
            dyn_phase.append(self._phases);
958
            dyn_time.append(0);
959
        
960
        step = time / steps;
961
        int_step = step / 10.0;
962
        
963
        for t in numpy.arange(step, time + step, step):
964
            # update states of oscillators
965
            self._phases = self._calculate_phases(solution, t, step, int_step);
966
            
967
            # update states of oscillators
968
            if (collect_dynamic == True):
969
                dyn_phase.append(self._phases);
970
                dyn_time.append(t);
971
        
972
        if (collect_dynamic != True):
973
            dyn_phase.append(self._phases);
974
            dyn_time.append(time);
975
                        
976
        output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time);
977
        return output_sync_dynamic;
978
979
980
    def _calculate_phases(self, solution, t, step, int_step):
981
        """!
982
        @brief Calculates new phases for oscillators in the network in line with current step.
983
        
984
        @param[in] solution (solve_type): Type solver of the differential equation.
985
        @param[in] t (double): Time of simulation.
986
        @param[in] step (double): Step of solution at the end of which states of oscillators should be calculated.
987
        @param[in] int_step (double): Step differentiation that is used for solving differential equation.
988
        
989
        @return (list) New states (phases) for oscillators.
990
        
991
        """
992
        
993
        next_phases = [0.0] * self._num_osc;    # new oscillator _phases
994
        
995
        for index in range (0, self._num_osc, 1):
996
            if (solution == solve_type.FAST):
997
                result = self._phases[index] + self._phase_kuramoto(self._phases[index], 0, index);
998
                next_phases[index] = self._phase_normalization(result);
999
                
1000
            elif (solution == solve_type.RK4):
1001
                result = odeint(self._phase_kuramoto, self._phases[index], numpy.arange(t - step, t, int_step), (index , ));
1002
                next_phases[index] = self._phase_normalization(result[len(result) - 1][0]);
1003
            
1004
            else:
1005
                raise NameError("Solver '" + solution + "' is not supported");
1006
        
1007
        return next_phases;
1008
1009
1010
    def _phase_normalization(self, teta):
1011
        """!
1012
        @brief Normalization of phase of oscillator that should be placed between [0; 2 * pi].
1013
        
1014
        @param[in] teta (double): phase of oscillator.
1015
        
1016
        @return (double) Normalized phase.
1017
        
1018
        """
1019
1020
        norm_teta = teta;
1021
        while (norm_teta > (2.0 * pi)) or (norm_teta < 0):
1022
            if (norm_teta > (2.0 * pi)):
1023
                norm_teta -= 2.0 * pi;
1024
            else:
1025
                norm_teta += 2.0 * pi;
1026
        
1027
        return norm_teta;
1028
1029
1030
    def get_neighbors(self, index):
1031
        """!
1032
        @brief Finds neighbors of the oscillator with specified index.
1033
        
1034
        @param[in] index (uint): index of oscillator for which neighbors should be found in the network.
1035
        
1036
        @return (list) Indexes of neighbors of the specified oscillator.
1037
        
1038
        """
1039
        
1040
        if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1041
            self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1042
            
1043
        return super().get_neighbors(index);
1044
1045
1046
    def has_connection(self, i, j):
1047
        """!
1048
        @brief Returns True if there is connection between i and j oscillators and False - if connection doesn't exist.
1049
        
1050
        @param[in] i (uint): index of an oscillator in the network.
1051
        @param[in] j (uint): index of an oscillator in the network.
1052
        
1053
        """
1054
        
1055
        if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1056
            self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1057
        
1058
        return super().has_connection(i, j);