diff --git a/hsyncnet/__init__.py b/hsyncnet/__init__.py
index b0b64f40..ebdfbe10 100644
--- a/hsyncnet/__init__.py
+++ b/hsyncnet/__init__.py
@@ -1,33 +1,70 @@
 from nnet.sync import *;
 
 from syncnet import syncnet;
-from support import average_neighbor_distance, read_sample, draw_clusters;
+from support import average_neighbor_distance, read_sample;
 
 class hsyncnet(syncnet):
     def __init__(self, source_data):
+        "Costructor of the oscillatory network hSync"
+        
+        "(in) source_data         - input data set defines structure of the network"
+        
         super().__init__(source_data);
     
     def process(self, number_clusters, order = 0.998, solution = solve_type.FAST, collect_dynamic = False):
+        "Perform clustering of input data set in line with input parameters"
+        
+        "(in) number_clusters     - number of clusters that should be allocated"
+        "(in) order               - level of local synchronization between oscillator that defines end of process of synchronization, range [0..1]"
+        "(in) solution            - type of solving differential equation: ode45, usual diff, etc."
+        "(in) collect_dynamic     - if True - return whole history of process synchronization otherwise - final state (when process of clustering is over)"
+        
+        "Returns dynamic of the network as tuple (time, oscillator_phases) that depends on collect_dynamic parameters"
+        
         number_neighbors = 3;
         current_number_clusters = numpy.Inf;
         
-        while(current_number_clusters > number_clusters):
-            radius = average_neighbor_distance(self._osc_loc, number_neighbors);
+        dyn_phase = [];
+        dyn_time = [];
+        
+        radius = average_neighbor_distance(self._osc_loc, number_neighbors);
+        
+        while(current_number_clusters > number_clusters):                
             self._create_connections(radius);
         
-            self.simulate_dynamic(order, solution, collect_dynamic);
-            clusters = self.get_clusters();
+            (t, dyn) = self.simulate_dynamic(order, solution, collect_dynamic);
+            if (collect_dynamic == True):
+                dyn_phase += dyn;
+                
+                if (len(dyn_time) > 0):
+                    point_time_last = dyn_time[len(dyn_time) - 1];
+                    dyn_time += [time_point + point_time_last for time_point in t];
+                else:
+                    dyn_time += t;
             
+            clusters = self.get_clusters(0.05);
+            
+            # Get current number of allocated clusters
             current_number_clusters = len(clusters);
+            
+            # Increase number of neighbors that should be used
             number_neighbors += 1;
+            
+            # Update connectivity radius and check if average function can be used anymore
+            if (number_neighbors >= len(self._osc_loc)):
+                radius = radius * 0.1 + radius;
+            else:
+                radius = average_neighbor_distance(self._osc_loc, number_neighbors);
         
+        
+        return (dyn_time, dyn_phase);
 
-sample = read_sample('../Samples/SampleSimple1.txt');
-network = hsyncnet(sample);
-
-network.process(2);
-
+# sample = read_sample('../Samples/SampleSimple1.txt');
+# network = hsyncnet(sample);
+#  
+# (dyn_time, dyn_phase) = network.process(2, collect_dynamic = True);
+#  
 # draw_dynamics(dyn_time, dyn_phase);
-
-clusters = network.get_clusters();
-draw_clusters(sample, clusters);
+#  
+# clusters = network.get_clusters();
+# draw_clusters(sample, clusters);
diff --git a/hsyncnet/examples.py b/hsyncnet/examples.py
new file mode 100644
index 00000000..9e1544f7
--- /dev/null
+++ b/hsyncnet/examples.py
@@ -0,0 +1,78 @@
+from nnet.sync import draw_dynamics;
+
+from support import read_sample, draw_clusters;
+from hsyncnet import hsyncnet;
+
+def template_clustering(file, number_clusters, arg_collect_dynamic = True, show_network_structure = False, arg_order = 0.999, arg_eps = 0.1):
+        sample = read_sample(file);
+        network = hsyncnet(sample);
+        
+        (time, dynamic) = network.process(number_clusters, order = arg_order, collect_dynamic = arg_collect_dynamic);
+        clusters = network.get_clusters();
+        
+        if (show_network_structure == True):
+            network.show_network();
+            
+        if (arg_collect_dynamic == True):
+            draw_dynamics(time, dynamic);
+        
+        draw_clusters(sample, clusters);
+        
+
+def cluster_sample1():
+    template_clustering('../Samples/SampleSimple1.txt', 2);
+    
+def cluster_sample2():
+    template_clustering('../Samples/SampleSimple2.txt', 3);
+    
+def cluster_sample3():
+    template_clustering('../Samples/SampleSimple3.txt', 4);
+    
+def cluster_simple4():
+    template_clustering('../Samples/SampleSimple4.txt', 5);
+    
+def cluster_simple5():
+    template_clustering('../Samples/SampleSimple5.txt', 4);
+    
+def cluster_elongate():
+    template_clustering('../Samples/SampleElongate.txt', 2, arg_collect_dynamic = False);
+
+def cluster_lsun():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleLsun.txt', 3, arg_collect_dynamic = False);
+    
+def cluster_hepta():
+    template_clustering('../Samples/SampleHepta.txt', 7, arg_collect_dynamic = False);
+    
+def cluster_tetra():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleTetra.txt', 4, arg_collect_dynamic = False);
+
+def cluster_target():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleTarget.txt', 6, arg_collect_dynamic = False);
+    
+def cluster_chainlink():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleChainlink.txt', 2, arg_collect_dynamic = False);
+    
+def cluster_wing_nut():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleWingNut.txt', 2, arg_collect_dynamic = False);
+    
+def cluster_two_diamonds():
+    "NOTE: Too slow"
+    template_clustering('../Samples/SampleTwoDiamonds.txt', 2, arg_collect_dynamic = False);    
+
+# cluster_sample1();
+# cluster_sample2();
+# cluster_sample3();
+# cluster_simple4();
+# cluster_elongate();
+# cluster_lsun();
+# cluster_hepta();
+# cluster_tetra();
+# cluster_target();
+# cluster_chainlink();
+# cluster_wing_nut();
+# cluster_two_diamonds();
diff --git a/hsyncnet/tests.py b/hsyncnet/tests.py
new file mode 100644
index 00000000..36ba1207
--- /dev/null
+++ b/hsyncnet/tests.py
@@ -0,0 +1,47 @@
+import unittest;
+
+from support import read_sample;
+from hsyncnet import hsyncnet;
+
+class Test(unittest.TestCase):
+    def templateClusteringResults(self, path, number_clusters, expected_length_clusters):
+        sample = read_sample(path);
+        network = hsyncnet(sample);
+        
+        (t, d) = network.process(number_clusters, order = 0.999, collect_dynamic = True);
+        clusters = network.get_clusters();
+        
+        assert sum([len(cluster) for cluster in clusters]) == sum(expected_length_clusters);
+        
+        if (sorted([len(cluster) for cluster in clusters]) != expected_length_clusters):
+#             print("Result: ", sorted([len(cluster) for cluster in clusters]), "Expect: ", expected_length_clusters);
+#             network.show_network();
+#             draw_dynamics(t, d);
+#             draw_clusters(sample, clusters);
+            
+            assert sorted([len(cluster) for cluster in clusters]) == expected_length_clusters;
+        
+    
+    def testClusteringSampleSimple1(self):
+        self.templateClusteringResults("../Samples/SampleSimple1.txt", 2, [5, 5]);
+        self.templateClusteringResults("../Samples/SampleSimple1.txt", 1, [10]);
+
+    def testClusteringSampleSimple2(self):
+        self.templateClusteringResults("../Samples/SampleSimple2.txt", 3, [5, 8, 10]);
+        self.templateClusteringResults("../Samples/SampleSimple2.txt", 1, [23]);
+
+    def testClusteringSampleSimple3(self):
+        self.templateClusteringResults("../Samples/SampleSimple3.txt", 4, [10, 10, 10, 30]);
+        self.templateClusteringResults("../Samples/SampleSimple3.txt", 1, [60]);
+        
+    def testClusteringSampleSimple4(self):
+        self.templateClusteringResults("../Samples/SampleSimple4.txt", 5, [15, 15, 15, 15, 15]);
+        self.templateClusteringResults("../Samples/SampleSimple4.txt", 1, [75]);
+        
+    def testClusteringSampleSimple5(self):
+        self.templateClusteringResults("../Samples/SampleSimple5.txt", 4, [15, 15, 15, 15]);
+        self.templateClusteringResults("../Samples/SampleSimple5.txt", 1, [60]);        
+
+
+if __name__ == "__main__":
+    unittest.main()
\ No newline at end of file
diff --git a/nnet/sync/__init__.py b/nnet/sync/__init__.py
index 5841d673..f3cc8cfb 100644
--- a/nnet/sync/__init__.py
+++ b/nnet/sync/__init__.py
@@ -9,7 +9,9 @@
 from scipy.integrate import odeint;
 from scipy.integrate import ode;
 from scipy.spatial import distance;
+
 from support import euclidean_distance;
+from hierarchical import hierarchical;
 
 
 class solve_type:
@@ -301,12 +303,9 @@ def simulate_dynamic(self, order = 0.998, solution = solve_type.FAST, collect_dy
         current_order = self.sync_local_order();
         
         # If requested input dynamics
-        dyn_phase = None;
-        dyn_time = None;
+        dyn_phase = [];
+        dyn_time = [];
         if (collect_dynamic == True):
-            dyn_phase = list();
-            dyn_time = list();
-            
             dyn_phase.append(self._phases);
             dyn_time.append(0);
         
diff --git a/syncnet/__init__.py b/syncnet/__init__.py
index cbe6c127..50192342 100644
--- a/syncnet/__init__.py
+++ b/syncnet/__init__.py
@@ -77,6 +77,12 @@ def get_neighbors(self, index):
     
     def phase_kuramoto(self, teta, t, argv):
         "Overrided method for calculation of oscillator phase"
+        
+        "(in) teta     - current value of phase"
+        "(in) t        - time (can be ignored)"
+        "(in) argv     - index of oscillator whose phase represented by argument teta"
+        
+        "Return new value of phase of oscillator with index 'argv'"
         index = argv;   # index of oscillator
         phase = 0;      # phase of a specified oscillator that will calculated in line with current env. states.
         
@@ -88,7 +94,11 @@ def phase_kuramoto(self, teta, t, argv):
 
 
     def get_clusters(self, eps = 0.1):
-        "Return clusters"
+        "Return list of clusters in line with state of ocillators (phases)."
+        
+        "(in) eps     - tolerance level that define maximal difference between phases of oscillators in one cluster"
+        
+        "Return list of clusters, for example [ [cluster1], [cluster2], ... ]"
         return self.allocate_sync_ensembles(eps);