GitHub

Harmonic Time Averages & Mesochronic Plots

Koopman-theoretic tools for revealing phase space structure through weighted time averages along trajectories.

Overview

Harmonic Time Averages (HTA) provide a Koopman-theoretic tool for analyzing the structure of phase space in discrete dynamical systems. By computing weighted time averages of observables along trajectories, HTAs reveal periodic orbits, resonance zones, and chaotic regions without requiring explicit integration of variational equations or Lyapunov exponent computation.

The core idea is simple: if a trajectory is periodic with a frequency that matches the test frequency, the harmonic average will be large in magnitude. If the trajectory is chaotic or periodic at an unrelated frequency, the average will tend toward zero as the number of iterations increases.

This implementation is based on the work of Mezic (2020) and Levnajic & Mezic (2014), which established the theoretical foundations of harmonic averages and mesochronic analysis for discrete-time systems.

Key references: I. Mezic, "Spectrum of the Koopman operator, spectral expansions in functional spaces, and state-space geometry," Journal of Nonlinear Science, 2020. Z. Levnajic and I. Mezic, "Ergodic theory and visualization," preprint, 2014.

HTA Definition

For an observable f, a map T, and a test frequency omega, the Harmonic Time Average is defined as:

HTA_omega(x) = lim_{N -> inf} (1/N) sum_{k=0}^{N-1} e^{2*pi*i*k*omega} f(T^k(x))

In practice, we compute a finite approximation by choosing a large but finite number of iterations N. The result is a complex number for each initial condition x. Two quantities carry dynamical information:

When the trajectory is chaotic or periodic at a frequency incommensurate with omega, the exponential weighting causes cancellation and the magnitude decays toward zero as N increases.

Observables

An observable is a scalar-valued function evaluated on the state vector at each iteration. The choice of observable affects the HTA result. The library provides several built-in observables:

Observable Definition Notes
Identity f(x) = x[0] Projects onto the first coordinate. Simple and general-purpose.
SinPi f(x) = sin(pi * x[0]) Smooth, bounded observable. Good for maps on [0, 1].
CosPi f(x) = cos(pi * x[0]) Complementary to SinPi. Even symmetry about x[0] = 0.
SinPiXY f(x) = sin(pi * x[0] * x[1]) Couples both coordinates. Useful for 2D maps.
Quadratic f(x) = x[0]^2 + x[1]^2 Radial observable. Measures distance from the origin.

You can also define custom observables by implementing the Observable trait (Rust), passing a callable (Python), or supplying a function (R).

Basic HTA Computation

The following example creates a standard map, picks an initial condition, computes the Harmonic Time Average at a chosen frequency, and inspects the magnitude and phase of the result.

use koopman_dmd::maps::StandardMap;
use koopman_dmd::harmonic::{harmonic_time_average, Observable};

fn main() {
    // Create a standard map with perturbation parameter K = 0.9
    let map = StandardMap::new(0.9);

    // Initial condition in phase space
    let x0 = [0.5, 0.25];

    // Test frequency (e.g., 1/3 for period-3 orbits)
    let omega = 1.0 / 3.0;

    // Compute HTA with 10000 iterations using the Identity observable
    let hta = harmonic_time_average(
        &map,
        &x0,
        omega,
        10_000,
        Observable::Identity,
    );

    // Inspect magnitude and phase
    println!("|HTA| = {:.6}", hta.norm());
    println!("arg(HTA) = {:.6}", hta.arg());
}
import koopman_dmd

# Create a standard map with perturbation parameter K = 0.9
map = koopman_dmd.StandardMap(k=0.9)

# Initial condition in phase space
x0 = [0.5, 0.25]

# Test frequency (e.g., 1/3 for period-3 orbits)
omega = 1.0 / 3.0

# Compute HTA with 10000 iterations using the Identity observable
hta = koopman_dmd.harmonic_time_average(
    map, x0, omega,
    n_iter=10_000,
    observable="identity",
)

# Inspect magnitude and phase
print(f"|HTA| = {abs(hta):.6f}")
print(f"arg(HTA) = {koopman_dmd.phase(hta):.6f}")
library(koopmandmd)

# Create a standard map with perturbation parameter K = 0.9
map <- standard_map(k = 0.9)

# Initial condition in phase space
x0 <- c(0.5, 0.25)

# Test frequency (e.g., 1/3 for period-3 orbits)
omega <- 1 / 3

# Compute HTA with 10000 iterations using the Identity observable
hta <- harmonic_time_average(
  map, x0, omega,
  n_iter = 10000,
  observable = "identity"
)

# Inspect magnitude and phase
cat(sprintf("|HTA| = %.6f\n", Mod(hta)))
cat(sprintf("arg(HTA) = %.6f\n", Arg(hta)))

Mesochronic Harmonic Plots

A mesochronic harmonic plot computes the HTA over a two-dimensional grid of initial conditions. Each pixel in the resulting image corresponds to a single initial condition; its color encodes the HTA magnitude (or phase). This produces a rich visualization of phase space structure: islands of stability appear as bright regions of high |HTA|, while chaotic seas appear dark.

The computation is embarrassingly parallel -- each initial condition is independent -- and the Rust implementation uses rayon for automatic work-stealing parallelism across all available CPU cores.

use koopman_dmd::maps::StandardMap;
use koopman_dmd::harmonic::{mesochronic_compute, Observable};

fn main() {
    let map = StandardMap::new(0.9);

    // Define the grid: 500x500 over [0, 2*pi) x [0, 2*pi)
    let resolution = 500;
    let x_range = (0.0, 2.0 * std::f64::consts::PI);
    let y_range = (0.0, 2.0 * std::f64::consts::PI);
    let omega = 1.0 / 3.0;
    let n_iter = 5_000;

    // Compute HTA over the grid (parallelized with rayon)
    let result = mesochronic_compute(
        &map,
        x_range,
        y_range,
        resolution,
        omega,
        n_iter,
        Observable::SinPi,
    );

    // result.hta_matrix: 500x500 Array2<f64> of |HTA| values
    // result.phase_matrix: 500x500 Array2<f64> of arg(HTA) values
    println!("HTA matrix shape: {:?}", result.hta_matrix.dim());
    println!("Phase matrix shape: {:?}", result.phase_matrix.dim());
}
import numpy as np
import koopman_dmd

map = koopman_dmd.StandardMap(k=0.9)

# Compute HTA over a 500x500 grid (parallelized internally via rayon)
result = koopman_dmd.mesochronic_compute(
    map,
    x_range=(0.0, 2.0 * np.pi),
    y_range=(0.0, 2.0 * np.pi),
    resolution=500,
    omega=1.0 / 3.0,
    n_iter=5_000,
    observable="sin_pi",
)

# result.hta_matrix: (500, 500) numpy array of |HTA| values
# result.phase_matrix: (500, 500) numpy array of arg(HTA) values
print("HTA matrix shape:", result.hta_matrix.shape)
print("Max |HTA|:", result.hta_matrix.max())
library(koopmandmd)

map <- standard_map(k = 0.9)

# Compute HTA over a 500x500 grid (parallelized internally via rayon)
result <- mesochronic_compute(
  map,
  x_range = c(0, 2 * pi),
  y_range = c(0, 2 * pi),
  resolution = 500,
  omega = 1 / 3,
  n_iter = 5000,
  observable = "sin_pi"
)

# result$hta_matrix: 500x500 matrix of |HTA| values
# result$phase_matrix: 500x500 matrix of arg(HTA) values
cat("HTA matrix dimensions:", dim(result$hta_matrix), "\n")
Interpretation: High |HTA| values indicate that the trajectory starting at that initial condition resonates with the test frequency omega. Resonance islands appear as bright connected regions. The chaotic sea, where orbits are ergodic and have continuous spectrum, appears uniformly dark.

Mesochronic Scatter Plots

The function mesochronic_scatter() computes HTA values for multiple observables simultaneously over the same grid of initial conditions. This is useful for producing scatter plots of |HTA_1| vs. |HTA_2|, where each point corresponds to one initial condition and can be colored by its position, phase, or classification.

Such scatter plots reveal clustering structure: initial conditions within the same resonance island tend to cluster together, while chaotic initial conditions cluster near the origin (both magnitudes small).

use koopman_dmd::maps::StandardMap;
use koopman_dmd::harmonic::{mesochronic_scatter, Observable};

fn main() {
    let map = StandardMap::new(0.9);

    // Compute HTA for two observables over a grid
    let result = mesochronic_scatter(
        &map,
        (0.0, 2.0 * std::f64::consts::PI),
        (0.0, 2.0 * std::f64::consts::PI),
        200,                     // grid resolution
        1.0 / 3.0,               // omega
        5_000,                   // iterations
        &[Observable::Identity, Observable::SinPi],
    );

    // result.hta_values: Vec<Vec<f64>> -- one entry per observable
    // result.initial_conditions: Vec<[f64; 2]> -- the grid points
    println!("Points: {}", result.initial_conditions.len());
    println!("Observables: {}", result.hta_values.len());
}
import numpy as np
import koopman_dmd

map = koopman_dmd.StandardMap(k=0.9)

# Compute HTA for two observables over a grid
result = koopman_dmd.mesochronic_scatter(
    map,
    x_range=(0.0, 2.0 * np.pi),
    y_range=(0.0, 2.0 * np.pi),
    resolution=200,
    omega=1.0 / 3.0,
    n_iter=5_000,
    observables=["identity", "sin_pi"],
)

# result.hta_values: list of numpy arrays, one per observable
# result.initial_conditions: (N, 2) array of grid points
print("Points:", result.initial_conditions.shape[0])
library(koopmandmd)

map <- standard_map(k = 0.9)

# Compute HTA for two observables over a grid
result <- mesochronic_scatter(
  map,
  x_range = c(0, 2 * pi),
  y_range = c(0, 2 * pi),
  resolution = 200,
  omega = 1 / 3,
  n_iter = 5000,
  observables = c("identity", "sin_pi")
)

# result$hta_values: list of numeric vectors, one per observable
# result$initial_conditions: Nx2 matrix of grid points
cat("Points:", nrow(result$initial_conditions), "\n")

Phase Space Classification

The function classify_phase_space() uses HTA values at multiple test frequencies to automatically classify each initial condition into one of three categories:

Category Criterion Interpretation
Resonating |HTA(omega)| > threshold The orbit is locked to the test frequency. It lies on or near a periodic orbit of the corresponding period.
NonResonating |HTA(omega)| < threshold but |HTA(omega')| is large for some other omega' The orbit is periodic, but at a frequency different from the test frequency. It belongs to a different island chain.
Chaotic |HTA| < threshold for all tested frequencies No persistent periodic structure detected. The orbit wanders ergodically through a chaotic region. 
use koopman_dmd::maps::StandardMap;
use koopman_dmd::harmonic::{classify_phase_space, Observable, Classification};

fn main() {
    let map = StandardMap::new(0.9);

    // Test frequencies: 1/2, 1/3, 1/4, 1/5
    let omegas = vec![0.5, 1.0/3.0, 0.25, 0.2];
    let threshold = 0.01;

    let result = classify_phase_space(
        &map,
        (0.0, 2.0 * std::f64::consts::PI),
        (0.0, 2.0 * std::f64::consts::PI),
        300,
        &omegas,
        10_000,
        threshold,
        Observable::Identity,
    );

    // Count classifications
    let n_resonating = result.iter()
        .filter(|c| matches!(c, Classification::Resonating(_)))
        .count();
    let n_chaotic = result.iter()
        .filter(|c| matches!(c, Classification::Chaotic))
        .count();

    println!("Resonating: {}", n_resonating);
    println!("Chaotic: {}", n_chaotic);
}
import numpy as np
import koopman_dmd

map = koopman_dmd.StandardMap(k=0.9)

# Test frequencies: 1/2, 1/3, 1/4, 1/5
omegas = [0.5, 1.0 / 3.0, 0.25, 0.2]
threshold = 0.01

result = koopman_dmd.classify_phase_space(
    map,
    x_range=(0.0, 2.0 * np.pi),
    y_range=(0.0, 2.0 * np.pi),
    resolution=300,
    omegas=omegas,
    n_iter=10_000,
    threshold=threshold,
    observable="identity",
)

# result.labels: array of strings ("resonating", "non_resonating", "chaotic")
# result.matched_omega: array of matched frequencies (NaN for chaotic)
print("Resonating:", sum(result.labels == "resonating"))
print("Chaotic:", sum(result.labels == "chaotic"))
library(koopmandmd)

map <- standard_map(k = 0.9)

# Test frequencies: 1/2, 1/3, 1/4, 1/5
omegas <- c(0.5, 1 / 3, 0.25, 0.2)
threshold <- 0.01

result <- classify_phase_space(
  map,
  x_range = c(0, 2 * pi),
  y_range = c(0, 2 * pi),
  resolution = 300,
  omegas = omegas,
  n_iter = 10000,
  threshold = threshold,
  observable = "identity"
)

# result$labels: character vector of classifications
# result$matched_omega: numeric vector of matched frequencies
table(result$labels)

HTA Convergence

The function hta_convergence() tracks how the HTA magnitude evolves as the number of iterations increases. This is useful for determining whether enough iterations were used to reliably distinguish resonant from non-resonant initial conditions.

For a resonating orbit, |HTA| converges to a positive constant. For a chaotic orbit, |HTA| decays roughly as 1/sqrt(N) due to the central limit theorem applied to the oscillating sum.

use koopman_dmd::maps::StandardMap;
use koopman_dmd::harmonic::{hta_convergence, Observable};

fn main() {
    let map = StandardMap::new(0.9);
    let x0 = [0.5, 0.25];
    let omega = 1.0 / 3.0;

    // Track |HTA| at checkpoints: 100, 500, 1000, 5000, 10000
    let checkpoints = vec![100, 500, 1_000, 5_000, 10_000];

    let curve = hta_convergence(
        &map,
        &x0,
        omega,
        &checkpoints,
        Observable::Identity,
    );

    // curve: Vec<(usize, f64)> -- (iterations, |HTA|)
    for (n, mag) in &curve {
        println!("N = {:>6} | |HTA| = {:.8}", n, mag);
    }
}
import koopman_dmd

map = koopman_dmd.StandardMap(k=0.9)
x0 = [0.5, 0.25]
omega = 1.0 / 3.0

# Track |HTA| at checkpoints
checkpoints = [100, 500, 1_000, 5_000, 10_000]

curve = koopman_dmd.hta_convergence(
    map, x0, omega,
    checkpoints=checkpoints,
    observable="identity",
)

# curve: list of (n_iter, |HTA|) pairs
for n, mag in curve:
    print(f"N = {n:>6} | |HTA| = {mag:.8f}")
library(koopmandmd)

map <- standard_map(k = 0.9)
x0 <- c(0.5, 0.25)
omega <- 1 / 3

# Track |HTA| at checkpoints
checkpoints <- c(100, 500, 1000, 5000, 10000)

curve <- hta_convergence(
  map, x0, omega,
  checkpoints = checkpoints,
  observable = "identity"
)

# curve: data.frame with columns n_iter and hta_magnitude
print(curve)

Mesochronic Sections

For higher-dimensional maps (such as the 4D Froeschle map), a full mesochronic plot over all dimensions is not practical. The function mesochronic_section() computes a two-dimensional slice through a higher-dimensional phase space: you fix two of the coordinates and vary the other two over a grid.

This is analogous to a Poincare section but computed via harmonic averages rather than intersection with a surface of section.

use koopman_dmd::maps::FroeschleMap;
use koopman_dmd::harmonic::{mesochronic_section, Observable};

fn main() {
    // 4D Froeschle map with coupling parameters
    let map = FroeschleMap::new(0.5, 0.3, 0.1);

    // Fix dimensions 2 and 3 at specific values
    let fixed_dims = vec![(2, 0.0), (3, 0.0)];

    // Vary dimensions 0 and 1 over the grid
    let vary_dims = (0, 1);
    let x_range = (0.0, 2.0 * std::f64::consts::PI);
    let y_range = (0.0, 2.0 * std::f64::consts::PI);

    let result = mesochronic_section(
        &map,
        vary_dims,
        x_range,
        y_range,
        &fixed_dims,
        300,          // resolution
        1.0 / 4.0,   // omega
        5_000,        // iterations
        Observable::Identity,
    );

    println!("Section shape: {:?}", result.hta_matrix.dim());
}
import numpy as np
import koopman_dmd

# 4D Froeschle map with coupling parameters
map = koopman_dmd.FroeschleMap(a=0.5, b=0.3, coupling=0.1)

# Fix dimensions 2 and 3, vary dimensions 0 and 1
result = koopman_dmd.mesochronic_section(
    map,
    vary_dims=(0, 1),
    x_range=(0.0, 2.0 * np.pi),
    y_range=(0.0, 2.0 * np.pi),
    fixed_dims={2: 0.0, 3: 0.0},
    resolution=300,
    omega=0.25,
    n_iter=5_000,
    observable="identity",
)

print("Section shape:", result.hta_matrix.shape)
library(koopmandmd)

# 4D Froeschle map with coupling parameters
map <- froeschle_map(a = 0.5, b = 0.3, coupling = 0.1)

# Fix dimensions 3 and 4, vary dimensions 1 and 2
result <- mesochronic_section(
  map,
  vary_dims = c(1, 2),
  x_range = c(0, 2 * pi),
  y_range = c(0, 2 * pi),
  fixed_dims = list("3" = 0.0, "4" = 0.0),
  resolution = 300,
  omega = 0.25,
  n_iter = 5000,
  observable = "identity"
)

cat("Section dimensions:", dim(result$hta_matrix), "\n")
Note: For the R bindings, dimension indices are 1-based (matching R convention), while the Rust and Python interfaces use 0-based indexing.