Van der Pol oscillator

examples/boundary_layer.py

@@ -35,7 +35,7 @@ def run_boundary_layer():
     K = 50 # Number of clusters
     L = 3 # Model order
 
-    # Create the Lorenz data
+    # Load the boundary layer data-set
     case_data = np.load('data/boundary_layer_l100_t120_a60.npz')
     data, dt = case_data['data'], case_data['dt']
     t = np.arange(data.shape[0]) * dt

examples/helper.py

@@ -320,7 +320,7 @@ def plot_cpd(x,x_hat):
     # Re-cluster original and CNM data with 10 clusters only for clarity
     from sklearn.cluster import KMeans
     K = 10
-    kmeans = KMeans(n_clusters=K,max_iter=300,n_init=10,n_jobs=-1)
+    kmeans = KMeans(n_clusters=K,max_iter=300,n_init=10)
     kmeans.fit(x)
     labels = kmeans.labels_
 
@@ -603,17 +603,17 @@ def lorenz(t,q,sigma,rho,beta):
     return data[points_to_remove:,:], dt
 
 def create_roessler_data():
-    """Create the Lorenz data"""
+    """Create the Roessler data"""
 
 
-    # Lorenz settings
+    # Roessler settings
     a = 0.1
     b = 0.1
     c = 14
     x0,y0,z0 = (1,1,1) # Initial conditions
     np.random.seed(0)
 
-    # Lorenz system
+    # Roessler system
     def Rossler(t,q,a,b,c):
         return [
                 -q[1] - q[2],
@@ -637,3 +637,54 @@ def Rossler(t,q,a,b,c):
 
     return data, t[1]-t[0]
 
+
+def create_vanderpol_data() -> tuple:
+    """Generates Van der Pol data.
+
+    Returns
+    -------
+    Tuple (data, dt). data is a ndarray containing the time-resolved data and 
+    dt is the float time-step value.
+    """
+
+    # Van der Pol settings
+    mu = 3.0                     # Van der Pol parameter "mu" 
+    initial_conditions = [1,1]   # Initial conditions x0, y0
+    np.random.seed(0)
+
+
+    def vanderpol(t: float, q: list, mu: float) -> list:
+        """Evaluates the Van der Pol oscillator equations.
+
+        Parameters
+        ----------
+        t:  time.
+        q:  list containing the variables [x, y].
+        mu: dissipation parameter.
+
+        Returns
+        -------
+        List containing the time derivatives [x', y'].
+        """
+        return  [
+                q[1],                                   # x' = y
+                mu * ( 1 - q[0]*q[0] ) * q[1] - q[0],   # y' = mu*(1-x^2)*y - x
+                ]
+
+    # Numerical settings
+    dt = 0.01             # Time step
+    N  = 50000            # Total number of steps
+    n_start = 3850        # Number of points to neglect at the beginning (to
+                          # remove the transient phase before the regular
+                          # oscillations)
+    N += n_start
+    t = np.arange(N) * dt # Time vector
+
+    # Integrate the ode
+    solution = solve_ivp(fun=lambda t, y: vanderpol(t, y, mu), 
+                         t_span = [0,t[-1]], 
+                         y0 = initial_conditions,
+                         t_eval=t)
+    data = solution.y[:,n_start:].T
+
+    return data, t[1]-t[0]

examples/kolmogorov.py

@@ -36,7 +36,7 @@ def run_kolmogorov():
     K = 200 # Number of clusters
     L = 24 # Model order
 
-    # Create the Lorenz data
+    # Load the Kolmogorov dataset
     case_data = np.load('data/kolmogorov.npz')
     data, dt = case_data['data'], case_data['dt']
     t = np.arange(data.shape[0]) * dt
@@ -45,7 +45,7 @@ def run_kolmogorov():
     # ----------
     cluster_config = {
             'data': data,
-            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10,n_jobs=-1),
+            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10),
             'dataset': 'kolmogorov',
             }
 

examples/lorenz.py

@@ -44,7 +44,7 @@ def run_lorenz():
     # ----------
     cluster_config = {
             'data': data,
-            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10,n_jobs=-1),
+            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10),
             'dataset': 'lorenz',
             }
 

examples/roessler.py

@@ -34,17 +34,17 @@ def run_roessler():
     # CNM parameters:
     # ---------------
     K = 100 # Number of clusters
-    L = 1 # Model order
+    L = 2 # Model order
 
-    # Create the Lorenz data
+    # Create the Roessler data
     data, dt = create_roessler_data()
     t = np.arange(data.shape[0]) * dt
 
     # Clustering
     # ----------
     cluster_config = {
             'data': data,
-            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10,n_jobs=-1),
+            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10),
             'dataset': 'roessler',
             }
 
@@ -85,7 +85,7 @@ def run_roessler():
     # time series
     time_range = (0,100)
     plot_label = ['x','y','z']
-    plot_time_series(t,data,t_hat,x_hat,time_range)
+    plot_time_series(t,data,t_hat,x_hat,time_range,plot_label)
 
     # cluster probability distribution
     plot_cpd(data,x_hat)

examples/vanderpol.py

@@ -0,0 +1,101 @@
+# -*- coding: utf-8 -*-
+#
+# Copyright (c) 2020 Daniel Fernex.
+# Copyright (c) 2020 Bernd R. Noack.
+# Copyright (c) 2020 Richard Semaan.
+#
+# This file is part of CNM 
+# (see https://github.com/fernexda/cnm).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+#
+
+import sys
+sys.path.insert(0,'../')
+from cnm import Clustering, TransitionProperties, Propagation
+import numpy as np
+
+def run_vanderpol():
+
+    from helper import create_vanderpol_data
+    from sklearn.cluster import KMeans
+
+    # CNM parameters:
+    # ---------------
+    K = 30 # Number of clusters
+    L = 3 # Model order
+
+    # Generate Van der Pol data:
+    data, dt = create_vanderpol_data()
+    t = np.arange(data.shape[0]) * dt
+
+    # Clustering
+    # ----------
+    cluster_config = {
+            'data': data,
+            'cluster_algo': KMeans(n_clusters=K,max_iter=300,n_init=10),
+            'dataset': 'vanderpol',
+            }
+
+    clustering = Clustering(**cluster_config)
+
+    # Transition properties
+    # ---------------------
+    transition_config = {
+            'clustering': clustering,
+            'dt': dt,
+            'K': K,
+            'L': L,
+            }
+
+    transition_properties = TransitionProperties(**transition_config)
+
+    # Propagation
+    # -----------
+    propagation_config = {
+            'transition_properties': transition_properties,
+            }
+
+    initial_centroid = 0 # Index of the centroid to start in
+    t_total = 950
+    dt_hat = dt      # To spline-interpolate the centroid-to-centroid trajectory
+
+    propagation = Propagation(**propagation_config)
+    t_hat, x_hat = propagation.run(t_total,initial_centroid,dt_hat)
+
+    # Plot the results
+    # ----------------
+    from helper import (plot_phase_space, plot_time_series,plot_cpd,
+                        plot_autocorrelation)
+
+    # phase space
+    n_dim = data.shape[1]
+    plot_phase_space(data,clustering.centroids,clustering.labels,n_dim=n_dim)
+
+    # time series
+    time_range = (45,60)
+    n_dim = 1
+    plot_label = ['x']
+    plot_time_series(t,data,t_hat,x_hat,time_range,plot_label,n_dim)
+
+    # cluster probability distribution
+    plot_cpd(data,x_hat)
+
+    # autocorrelation function
+    time_blocks = 40
+    time_range = (-0.5,time_blocks)
+    plot_autocorrelation(t,data,t_hat,x_hat,time_blocks,time_range)
+
+if __name__== '__main__':
+    run_vanderpol()