Completed
Push — 0.7.dev ( 8e244d...3d4397 )
by Andrei
01:10
created

sync_dynamic.calculate_sync_order()   B

Complexity

Conditions 2

Size

Total Lines 26

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 2
dl 0
loc 26
rs 8.8571
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
        sequence_order = [];
368
        for index in range(start_iteration, stop_iteration):
369
            sequence_order.append(order_estimator.calculate_sync_order(self.output[index]));
370
        
371
        return sequence_order;
372
373
374
    def calculate_local_order_parameter(self, oscillatory_network, start_iteration = None, stop_iteration = None):
375
        """!
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...
376
        @brief Calculates local order parameter.
377
        @details Local order parameter or so-called level of local or partial synchronization is calculated by following expression:
378
        
379
        \f[
380 View Code Duplication
        r_{c}=\left | \sum_{i=0}^{N} \frac{1}{N_{i}} \sum_{j=0}e^{ \theta_{j} - \theta_{i} } \right |;
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
381
        \f]
382
        
383
        where N - total amount of oscillators in the network and \f$N_{i}\f$ - amount of neighbors of oscillator with index \f$i\f$.
384
        
385
        @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
386
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the last iteration is used.
387
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then 'start_iteration' + 1 is used.
388
        
389
        @return (list) List of levels of local (partial) synchronization (local order parameter evolution).
390
        
391
        """
392
393
        (start_iteration, stop_iteration) = self.__get_start_stop_iterations(start_iteration, stop_iteration);
394
        
395
        sequence_local_order = [];
396
        for index in range(start_iteration, stop_iteration):
397
            sequence_local_order.append(order_estimator.calculate_local_sync_order(self.output[index], oscillatory_network));
398
        
399
        return sequence_local_order;
400
401
402
    def __get_start_stop_iterations(self, start_iteration, stop_iteration):
403
        """!
404
        @brief Aplly rules for start_iteration and stop_iteration parameters.
405
406
        @param[in] start_iteration (uint): The first iteration that is used for calculation.
407
        @param[in] stop_iteration (uint): The last iteration that is used for calculation.
408
        
409
        @return (tuple) New the first iteration and the last.
410
        
411
        """
412
        if (start_iteration is None):
413
            start_iteration = len(self._dynamic) - 1;
414 View Code Duplication
        
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
415
        if (stop_iteration is None):
416
            stop_iteration = start_iteration + 1;
417
        
418
        return (start_iteration, stop_iteration);
419
420
421
422
class sync_visualizer:
423
    """!
424
    @brief Visualizer of output dynamic of sync network (Sync).
425
    
426
    """
427
428
    @staticmethod
429
    def show_output_dynamic(sync_output_dynamic):
430
        """!
431
        @brief Shows output dynamic (output of each oscillator) during simulation.
432
        
433
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
434
        
435
        @see show_output_dynamics
436
        
437
        """
438
        
439
        draw_dynamics(sync_output_dynamic.time, sync_output_dynamic.output, x_title = "t", y_title = "phase", y_lim = [0, 2 * 3.14]);
440
    
441
    
442
    @staticmethod
443
    def show_output_dynamics(sync_output_dynamics):
444
        """!
445
        @brief Shows several output dynamics (output of each oscillator) during simulation.
446
        @details Each dynamic is presented on separate plot.
447
        
448
        @param[in] sync_output_dynamics (list): list of output dynamics 'sync_dynamic' of the Sync network.
449
        
450
        @see show_output_dynamic
451
        
452
        """
453
        
454
        draw_dynamics_set(sync_output_dynamics, "t", "phase", None, [0, 2 * 3.14], False, False);
455
    
456
    
457
    @staticmethod
458
    def show_correlation_matrix(sync_output_dynamic, iteration = None):
459
        """!
460
        @brief Shows correlation matrix between oscillators at the specified iteration.
461
        
462
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
463
        @param[in] iteration (uint): Number of interation of simulation for which correlation matrix should be allocated.
464
                                      If iternation number is not specified, the last step of simulation is used for the matrix allocation.
465
        
466
        """
467
        
468
        _ = plt.figure();
469
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(iteration);
470
        
471
        plt.imshow(correlation_matrix, cmap = plt.get_cmap('cool'), interpolation='kaiser', vmin = 0.0, vmax = 1.0); 
472
        plt.show();
473
474
475
    @staticmethod
476
    def show_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, iteration = None):
477
        """!
478
        @brief Shows 2D matrix of phase values of oscillators at the specified iteration.
479
        @details User should ensure correct matrix sizes in line with following expression grid_width x grid_height that should be equal to 
480
                  amount of oscillators otherwise exception is thrown. If grid_width or grid_height are not specified than phase matrix size 
481
                  will by calculated automatically by square root.
482
        
483
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose phase matrix should be shown.
484
        @param[in] grid_width (uint): Width of the phase matrix.
485
        @param[in] grid_height (uint): Height of the phase matrix.
486
        @param[in] iteration (uint): Number of iteration of simulation for which correlation matrix should be allocated.
487
                    If iternation number is not specified, the last step of simulation is used for the matrix allocation.
488
        
489
        """
490
        
491
        _ = plt.figure();
492
        phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, iteration);
493
        
494
        plt.imshow(phase_matrix, cmap = plt.get_cmap('jet'), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi); 
495
        plt.show();
496
497
498
    @staticmethod
499
    def show_order_parameter(sync_output_dynamic, start_iteration = None, stop_iteration = None):
500
        """!
501
        @brief Shows evolution of order parameter (level of global synchronization in the network).
502
        
503
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
504
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
505
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
506
        
507
        """
508
        
509
        (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
510
        
511
        order_parameter = sync_output_dynamic.calculate_order_parameter(start_iteration, stop_iteration);
512
        axis = plt.subplot(111);
513
        plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
514
        set_ax_param(axis, "t", "R (order parameter)", None, [0.0, 1.05]);
515
        
516
        plt.show();
517
518
519
    @staticmethod
520
    def show_local_order_parameter(sync_output_dynamic, oscillatory_network, start_iteration = None, stop_iteration = None):
521
        """!
522
        @brief Shows evolution of local order parameter (level of local synchronization in the network).
523
        
524
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network whose evolution of global synchronization should be visualized.
525
        @param[in] oscillatory_network (sync): Sync oscillatory network whose structure of connections is required for calculation.
526
        @param[in] start_iteration (uint): The first iteration that is used for calculation, if 'None' then the first is used
527
        @param[in] stop_iteration (uint): The last iteration that is used for calculation, if 'None' then the last is used.
528
        
529
        """
530
        (start_iteration, stop_iteration) = sync_visualizer.__get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration);
531
        
532
        order_parameter = sync_output_dynamic.calculate_local_order_parameter(oscillatory_network, start_iteration, stop_iteration);
533
        axis = plt.subplot(111);
534
        plt.plot(sync_output_dynamic.time[start_iteration:stop_iteration], order_parameter, 'b-', linewidth = 2.0);
535
        set_ax_param(axis, "t", "R (local order parameter)", None, [0.0, 1.05]);
536
        
537
        plt.show();
538
539
540
    @staticmethod
541
    def animate_output_dynamic(sync_output_dynamic, animation_velocity = 75, save_movie = None):
542
        """!
543
        @brief Shows animation of output dynamic (output of each oscillator) during simulation on a circle from [0; 2pi].
544
        
545
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
546
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
547
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
548
        
549
        """
550
        
551
        figure = plt.figure();
552
        
553
        dynamic = sync_output_dynamic.output[0];
554
        artist, = plt.polar(dynamic, [1.0] * len(dynamic), 'o', color = 'blue');
555
        
556
        def init_frame():
557
            return [ artist ];
558
        
559
        def frame_generation(index_dynamic):
560
            dynamic = sync_output_dynamic.output[index_dynamic];
561
            artist.set_data(dynamic, [1.0] * len(dynamic));
562
            
563
            return [ artist ];
564
        
565
        phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = animation_velocity, init_func = init_frame, repeat_delay = 5000);
566
567
        if (save_movie is not None):
568
            phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
569
        else:
570
            plt.show();
571
572
573
    @staticmethod
574
    def animate_correlation_matrix(sync_output_dynamic, animation_velocity = 75, colormap = 'cool', save_movie = None):
575
        """!
576
        @brief Shows animation of correlation matrix between oscillators during simulation.
577
        
578
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
579
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
580
        @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
581
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
582
        
583
        """
584
        
585
        figure = plt.figure();
586
        
587
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0);
588
        artist = plt.imshow(correlation_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', hold = True, vmin = 0.0, vmax = 1.0);
589
        
590
        def init_frame(): 
591
            return [ artist ];
592
        
593
        def frame_generation(index_dynamic):
594
            correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic);
595
            artist.set_data(correlation_matrix);
596
            
597
            return [ artist ];
598
599
        correlation_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000, blit = True);
600
        
601
        if (save_movie is not None):
602
            correlation_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
603
        else:
604
            plt.show();
605
606
607
    @staticmethod
608
    def animate_phase_matrix(sync_output_dynamic, grid_width = None, grid_height = None, animation_velocity = 75, colormap = 'jet', save_movie = None):
609
        """!
610
        @brief Shows animation of phase matrix between oscillators during simulation on 2D stage.
611
        @details If grid_width or grid_height are not specified than phase matrix size will by calculated automatically by square root.
612
        
613
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
614
        @param[in] grid_width (uint): Width of the phase matrix.
615
        @param[in] grid_height (uint): Height of the phase matrix.
616
        @param[in] animation_velocity (uint): Interval between frames in milliseconds.
617
        @param[in] colormap (string): Name of colormap that is used by matplotlib ('gray', 'pink', 'cool', spring', etc.).
618
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
619
        
620
        """
621
        
622
        figure = plt.figure();
623
        
624
        def init_frame(): 
625
            return frame_generation(0);
626
        
627
        def frame_generation(index_dynamic):
628
            figure.clf();
629
            axis = figure.add_subplot(111);
630
            
631
            phase_matrix = sync_output_dynamic.allocate_phase_matrix(grid_width, grid_height, index_dynamic);
632
            axis.imshow(phase_matrix, cmap = plt.get_cmap(colormap), interpolation='kaiser', vmin = 0.0, vmax = 2.0 * math.pi);
633
            artist = figure.gca();
634
            
635
            return [ artist ];
636
637
        phase_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), init_func = init_frame, interval = animation_velocity , repeat_delay = 1000);
638
        
639
        if (save_movie is not None):
640
            phase_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
641
        else:
642
            plt.show();
643
644
645
    @staticmethod
646
    def __get_start_stop_iterations(sync_output_dynamic, start_iteration, stop_iteration):
647
        """!
648
        @brief Apply rule of preparation for start iteration and stop iteration values.
649
        
650
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
651
        @param[in] start_iteration (uint): The first iteration that is used for calculation.
652
        @param[in] stop_iteration (uint): The last iteration that is used for calculation.
653
        
654
        @return (tuple) New values of start and stop iterations.
655
        
656
        """
657
        if (start_iteration is None):
658
            start_iteration = 0;
659
        
660
        if (stop_iteration is None):
661
            stop_iteration = len(sync_output_dynamic);
662
        
663
        return (start_iteration, stop_iteration);
664
665
666
    @staticmethod
667
    def animate(sync_output_dynamic, title = None, save_movie = None):
668
        """!
669
        @brief Shows animation of phase coordinates and animation of correlation matrix together for the Sync dynamic output on the same figure.
670
        
671
        @param[in] sync_output_dynamic (sync_dynamic): Output dynamic of the Sync network.
672
        @param[in] title (string): Title of the animation that is displayed on a figure if it is specified.
673
        @param[in] save_movie (string): If it is specified then animation will be stored to file that is specified in this parameter.
674
        
675
        """
676
        
677
        dynamic = sync_output_dynamic.output[0];
678
        correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(0);
679
        
680
        figure = plt.figure(1);
681
        if (title is not None):
682
            figure.suptitle(title, fontsize = 26, fontweight = 'bold')
683
        
684
        ax1 = figure.add_subplot(121, projection='polar');
685
        ax2 = figure.add_subplot(122);
686
        
687
        artist1, = ax1.plot(dynamic, [1.0] * len(dynamic), marker = 'o', color = 'blue', ls = '');
688
        artist2 = ax2.imshow(correlation_matrix, cmap = plt.get_cmap('Accent'), interpolation='kaiser');
689
        
690
        def init_frame():
691
            return [ artist1, artist2 ];
692
693
        def frame_generation(index_dynamic):
694
            dynamic = sync_output_dynamic.output[index_dynamic];
695
            artist1.set_data(dynamic, [1.0] * len(dynamic));
696
            
697
            correlation_matrix = sync_output_dynamic.allocate_correlation_matrix(index_dynamic);
698
            artist2.set_data(correlation_matrix);
699
            
700
            return [ artist1, artist2 ];
701
        
702
        dynamic_animation = animation.FuncAnimation(figure, frame_generation, len(sync_output_dynamic), interval = 75, init_func = init_frame, repeat_delay = 5000);
703
        
704
        if (save_movie is not None):
705
            dynamic_animation.save(save_movie, writer = 'ffmpeg', fps = 15, bitrate = 1500);
706
        else:
707
            plt.show();
708
709
710
711
class sync_network(network):
712
    """!
713
    @brief Model of oscillatory network that is based on the Kuramoto model of synchronization.
714
    
715
    """
716
717
    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):
718
        """!
719
        @brief Constructor of oscillatory network is based on Kuramoto model.
720
        
721
        @param[in] num_osc (uint): Number of oscillators in the network.
722
        @param[in] weight (double): Coupling strength of the links between oscillators.
723
        @param[in] frequency (double): Multiplier of internal frequency of the oscillators.
724
        @param[in] type_conn (conn_type): Type of connection between oscillators in the network (all-to-all, grid, bidirectional list, etc.).
725
        @param[in] representation (conn_represent): Internal representation of connection in the network: matrix or list.
726
        @param[in] initial_phases (initial_type): Type of initialization of initial phases of oscillators (random, uniformly distributed, etc.).
727
        @param[in] ccore (bool): If True simulation is performed by CCORE library (C++ implementation of pyclustering).
728
        
729
        """
730
        
731
        self._ccore_network_pointer = None;      # Pointer to CCORE Sync implementation of the network.
732
        
733
        if (ccore is True):
734
            self._ccore_network_pointer = wrapper.sync_create_network(num_osc, weight, frequency, type_conn, initial_phases);
735
            self._num_osc = num_osc;
736
            self._conn_represent = conn_represent.MATRIX;
737
        
738
        else:   
739
            super().__init__(num_osc, type_conn, representation);
740
            
741
            self._weight = weight;
742
            
743
            self._phases = list();
744
            self._freq = list();
745
            
746
            random.seed();
747
            for index in range(0, num_osc, 1):
748
                if (initial_phases == initial_type.RANDOM_GAUSSIAN):
749
                    self._phases.append(random.random() * 2 * pi);
750
                
751
                elif (initial_phases == initial_type.EQUIPARTITION):
752
                    self._phases.append( pi / num_osc * index);
753
                
754
                self._freq.append(random.random() * frequency);
755
756
757
    def __del__(self):
758
        """!
759
        @brief Destructor of oscillatory network is based on Kuramoto model.
760
        
761
        """
762
        
763
        if (self._ccore_network_pointer is not None):
764
            wrapper.sync_destroy_network(self._ccore_network_pointer);
765
            self._ccore_network_pointer = None;
766
767
768
    def sync_order(self):
769
        """!
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...
770
        @brief Calculates current level of global synchorization (order parameter) in the network.
771
        @details This parameter is tend 1.0 when the oscillatory network close to global synchronization and it tend to 0.0 when 
772
                  desynchronization is observed in the network. Order parameter is calculated using following equation:
773
                  
774
                  \f[
775
                  r_{c}=\frac{1}{Ne^{i\varphi }}\sum_{j=0}^{N}e^{i\theta_{j}};
776
                  \f]
777
                  
778
                  where \f$\varphi\f$ is a average phase coordinate in the network, \f$N\f$ is an amount of oscillators in the network.
779
        
780
        Example:
781
        @code
782
            oscillatory_network = sync(16, type_conn = conn_type.ALL_TO_ALL);
783
            output_dynamic = oscillatory_network.simulate_static(100, 10);
784
            
785
            if (oscillatory_network.sync_order() < 0.9): print("Global synchronization is not reached yet.");
786
            else: print("Global synchronization is reached.");
787
        @endcode
788
        
789
        @return (double) Level of global synchronization (order parameter).
790
        
791
        @see sync_local_order()
792
        
793
        """
794
        
795
        if (self._ccore_network_pointer is not None):
796
            return wrapper.sync_order(self._ccore_network_pointer);
797
        
798
        return order_estimator.calculate_sync_order(self._phases);
799
800
801
    def sync_local_order(self):
802
        """!
803
        @brief Calculates current level of local (partial) synchronization in the network.
804
        
805
        @return (double) Level of local (partial) synchronization.
806
        
807
        @see sync_order()
808
        
809
        """
810
        
811
        if (self._ccore_network_pointer is not None):
812
            return wrapper.sync_local_order(self._ccore_network_pointer);
813
        
814
        return order_estimator.calculate_local_sync_order(self._phases, self);
815
816
817
    def _phase_kuramoto(self, teta, t, argv):
0 ignored issues
show
Unused Code introduced by
The argument t seems to be unused.
Loading history...
818
        """!
819
        @brief Returns result of phase calculation for specified oscillator in the network.
820
        
821
        @param[in] teta (double): Phase of the oscillator that is differentiated.
822
        @param[in] t (double): Current time of simulation.
823
        @param[in] argv (tuple): Index of the oscillator in the list.
824
        
825
        @return (double) New phase for specified oscillator (don't assign here).
826
        
827
        """
828
        
829
        index = argv;
830
        phase = 0;
831
        for k in range(0, self._num_osc):
832
            if (self.has_connection(index, k) == True):
833
                phase += math.sin(self._phases[k] - teta);
834
            
835
        return ( self._freq[index] + (phase * self._weight / self._num_osc) );
836
837
838
    def simulate(self, steps, time, solution = solve_type.FAST, collect_dynamic = True):
839
        """!
840
        @brief Performs static simulation of Sync oscillatory network.
841
        
842
        @param[in] steps (uint): Number steps of simulations during simulation.
843
        @param[in] time (double): Time of simulation.
844
        @param[in] solution (solve_type): Type of solution (solving).
845
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
846
        
847
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
848
                otherwise returns only last values (last step of simulation) of dynamic.
849
        
850
        @see simulate_dynamic()
851
        @see simulate_static()
852
        
853
        """
854
        
855
        return self.simulate_static(steps, time, solution, collect_dynamic);
856
857
858
    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):
859
        """!
860
        @brief Performs dynamic simulation of the network until stop condition is not reached. Stop condition is defined by input argument 'order'.
861
        
862
        @param[in] order (double): Order of process synchronization, distributed 0..1.
863
        @param[in] solution (solve_type): Type of solution.
864
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
865
        @param[in] step (double): Time step of one iteration of simulation.
866
        @param[in] int_step (double): Integration step, should be less than step.
867
        @param[in] threshold_changes (double): Additional stop condition that helps prevent infinite simulation, defines limit of changes of oscillators between current and previous steps.
868
        
869
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
870
                otherwise returns only last values (last step of simulation) of dynamic.
871
        
872
        @see simulate()
873
        @see simulate_static()
874
        
875
        """
876
        
877
        if (self._ccore_network_pointer is not None):
878
            ccore_instance_dynamic = wrapper.sync_simulate_dynamic(self._ccore_network_pointer, order, solution, collect_dynamic, step, int_step, threshold_changes);
879
            return sync_dynamic(None, None, ccore_instance_dynamic);
880
        
881
        # For statistics and integration
882
        time_counter = 0;
883
        
884
        # Prevent infinite loop. It's possible when required state cannot be reached.
885
        previous_order = 0;
886
        current_order = self.sync_local_order();
887
        
888
        # If requested input dynamics
889
        dyn_phase = [];
890
        dyn_time = [];
891
        if (collect_dynamic == True):
892
            dyn_phase.append(self._phases);
893
            dyn_time.append(0);
894
        
895
        # Execute until sync state will be reached
896
        while (current_order < order):
897
            # update states of oscillators
898
            self._phases = self._calculate_phases(solution, time_counter, step, int_step);
899
            
900
            # update time
901
            time_counter += step;
902
            
903
            # if requested input dynamic
904
            if (collect_dynamic == True):
905
                dyn_phase.append(self._phases);
906
                dyn_time.append(time_counter);
907
                
908
            # update orders
909
            previous_order = current_order;
910
            current_order = self.sync_local_order();
911
            
912
            # hang prevention
913
            if (abs(current_order - previous_order) < threshold_changes):
914
                # print("Warning: sync_network::simulate_dynamic - simulation is aborted due to low level of convergence rate (order = " + str(current_order) + ").");
915
                break;
916
            
917
        if (collect_dynamic != True):
918
            dyn_phase.append(self._phases);
919
            dyn_time.append(time_counter);
920
921
        output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time, None);
922
        return output_sync_dynamic;
923
924
925
    def simulate_static(self, steps, time, solution = solve_type.FAST, collect_dynamic = False):
926
        """!
927
        @brief Performs static simulation of oscillatory network.
928
        
929
        @param[in] steps (uint): Number steps of simulations during simulation.
930
        @param[in] time (double): Time of simulation.
931
        @param[in] solution (solve_type): Type of solution.
932
        @param[in] collect_dynamic (bool): If True - returns whole dynamic of oscillatory network, otherwise returns only last values of dynamics.
933
        
934
        @return (list) Dynamic of oscillatory network. If argument 'collect_dynamic' = True, than return dynamic for the whole simulation time,
935
                otherwise returns only last values (last step of simulation) of dynamic.
936
        
937
        @see simulate()
938
        @see simulate_dynamic()
939
        
940
        """
941
        
942
        if (self._ccore_network_pointer is not None):
943
            ccore_instance_dynamic = wrapper.sync_simulate_static(self._ccore_network_pointer, steps, time, solution, collect_dynamic);
944
            return sync_dynamic(None, None, ccore_instance_dynamic);
945
        
946
        dyn_phase = [];
947
        dyn_time = [];
948
        
949
        if (collect_dynamic == True):
950
            dyn_phase.append(self._phases);
951
            dyn_time.append(0);
952
        
953
        step = time / steps;
954
        int_step = step / 10.0;
955
        
956
        for t in numpy.arange(step, time + step, step):
957
            # update states of oscillators
958
            self._phases = self._calculate_phases(solution, t, step, int_step);
959
            
960
            # update states of oscillators
961
            if (collect_dynamic == True):
962
                dyn_phase.append(self._phases);
963
                dyn_time.append(t);
964
        
965
        if (collect_dynamic != True):
966
            dyn_phase.append(self._phases);
967
            dyn_time.append(time);
968
                        
969
        output_sync_dynamic = sync_dynamic(dyn_phase, dyn_time);
970
        return output_sync_dynamic;
971
972
973
    def _calculate_phases(self, solution, t, step, int_step):
974
        """!
975
        @brief Calculates new phases for oscillators in the network in line with current step.
976
        
977
        @param[in] solution (solve_type): Type solver of the differential equation.
978
        @param[in] t (double): Time of simulation.
979
        @param[in] step (double): Step of solution at the end of which states of oscillators should be calculated.
980
        @param[in] int_step (double): Step differentiation that is used for solving differential equation.
981
        
982
        @return (list) New states (phases) for oscillators.
983
        
984
        """
985
        
986
        next_phases = [0.0] * self._num_osc;    # new oscillator _phases
987
        
988
        for index in range (0, self._num_osc, 1):
989
            if (solution == solve_type.FAST):
990
                result = self._phases[index] + self._phase_kuramoto(self._phases[index], 0, index);
991
                next_phases[index] = self._phase_normalization(result);
992
                
993
            elif (solution == solve_type.RK4):
994
                result = odeint(self._phase_kuramoto, self._phases[index], numpy.arange(t - step, t, int_step), (index , ));
995
                next_phases[index] = self._phase_normalization(result[len(result) - 1][0]);
996
            
997
            else:
998
                raise NameError("Solver '" + solution + "' is not supported");
999
        
1000
        return next_phases;
1001
1002
1003
    def _phase_normalization(self, teta):
1004
        """!
1005
        @brief Normalization of phase of oscillator that should be placed between [0; 2 * pi].
1006
        
1007
        @param[in] teta (double): phase of oscillator.
1008
        
1009
        @return (double) Normalized phase.
1010
        
1011
        """
1012
1013
        norm_teta = teta;
1014
        while (norm_teta > (2.0 * pi)) or (norm_teta < 0):
1015
            if (norm_teta > (2.0 * pi)):
1016
                norm_teta -= 2.0 * pi;
1017
            else:
1018
                norm_teta += 2.0 * pi;
1019
        
1020
        return norm_teta;
1021
1022
1023
    def get_neighbors(self, index):
1024
        """!
1025
        @brief Finds neighbors of the oscillator with specified index.
1026
        
1027
        @param[in] index (uint): index of oscillator for which neighbors should be found in the network.
1028
        
1029
        @return (list) Indexes of neighbors of the specified oscillator.
1030
        
1031
        """
1032
        
1033
        if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1034
            self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1035
            
1036
        return super().get_neighbors(index);
1037
1038
1039
    def has_connection(self, i, j):
1040
        """!
1041
        @brief Returns True if there is connection between i and j oscillators and False - if connection doesn't exist.
1042
        
1043
        @param[in] i (uint): index of an oscillator in the network.
1044
        @param[in] j (uint): index of an oscillator in the network.
1045
        
1046
        """
1047
        
1048
        if ( (self._ccore_network_pointer is not None) and (self._osc_conn is None) ):
1049
            self._osc_conn = wrapper.sync_connectivity_matrix(self._ccore_network_pointer);
1050
        
1051
        return super().has_connection(i, j);