From b112531eedc37240101ceddd663b45e18a8c815f Mon Sep 17 00:00:00 2001 From: Paolo Olivucci Date: Wed, 24 Nov 2021 13:40:53 +0000 Subject: [PATCH 1/2] added Van der Pol oscillator, fixed typos --- examples/boundary_layer.py | 2 +- examples/helper.py | 59 +++++++++++++++++++-- examples/kolmogorov.py | 4 +- examples/lorenz.py | 2 +- examples/roessler.py | 8 +-- examples/vanderpol.py | 102 +++++++++++++++++++++++++++++++++++++ 6 files changed, 165 insertions(+), 12 deletions(-) create mode 100644 examples/vanderpol.py diff --git a/examples/boundary_layer.py b/examples/boundary_layer.py index ddb9d48..e50fc9a 100644 --- a/examples/boundary_layer.py +++ b/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 diff --git a/examples/helper.py b/examples/helper.py index 99014fa..a02c1a0 100644 --- a/examples/helper.py +++ b/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] \ No newline at end of file diff --git a/examples/kolmogorov.py b/examples/kolmogorov.py index 3bc3251..ea89fe6 100644 --- a/examples/kolmogorov.py +++ b/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', } diff --git a/examples/lorenz.py b/examples/lorenz.py index 51d16c6..aa2a176 100644 --- a/examples/lorenz.py +++ b/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', } diff --git a/examples/roessler.py b/examples/roessler.py index 2e93301..035ea53 100644 --- a/examples/roessler.py +++ b/examples/roessler.py @@ -34,9 +34,9 @@ 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 @@ -44,7 +44,7 @@ def run_roessler(): # ---------- 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) diff --git a/examples/vanderpol.py b/examples/vanderpol.py new file mode 100644 index 0000000..46a7833 --- /dev/null +++ b/examples/vanderpol.py @@ -0,0 +1,102 @@ +# -*- coding: utf-8 -*- +# +# Copyright (c) 2021 Paolo Olivucci. +# 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 . +# + +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() From 27e608a3b074ddeca5776708025e233635957492 Mon Sep 17 00:00:00 2001 From: Paolo Olivucci Date: Fri, 26 Nov 2021 15:16:00 +0000 Subject: [PATCH 2/2] Fixed typos --- examples/vanderpol.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/vanderpol.py b/examples/vanderpol.py index 46a7833..832186b 100644 --- a/examples/vanderpol.py +++ b/examples/vanderpol.py @@ -1,6 +1,5 @@ # -*- coding: utf-8 -*- # -# Copyright (c) 2021 Paolo Olivucci. # Copyright (c) 2020 Daniel Fernex. # Copyright (c) 2020 Bernd R. Noack. # Copyright (c) 2020 Richard Semaan.