From ffc1010f5c2dc4c7559677bd8f5e827beea90df7 Mon Sep 17 00:00:00 2001
From: Claude Code <noreply@anthropic.com>
Date: Wed, 7 Jan 2026 19:41:49 -0500
Subject: [PATCH] Add Smooth Orthogonal Decomposition (SOD) implementation
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Implements the Smooth Orthogonal Decomposition (SOD) algorithm for extracting
linear normal modes and natural frequencies from time-series data, as proposed
by Chelidze and Zhou (2006).

SOD extends POD by considering temporal smoothness in addition to spatial variance.
It uses the covariance matrices of displacement and velocity responses to form a
generalized eigenvalue problem (Σ_xx Ψ = λ Σ_vv Ψ), which yields smooth orthogonal
modes and values without requiring knowledge of the system's mass matrix.

Key features:
- SOD type following the same pattern as POD
- Automatic velocity computation from displacement data using finite differences
- Support for user-provided velocity data
- Energy-based truncation for automatic mode selection
- Natural frequency estimation from smooth orthogonal values (SOVs)
- Support for both matrix and vector-of-vectors input formats

References:
- Chelidze D, Zhou W. Smooth orthogonal decomposition-based vibration mode
  identification. J Sound Vib. 2006; 292(3-5): 461-473.

Closes #9

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
---
 src/DataReduction/SOD.jl   | 243 +++++++++++++++++++++++++++++++++++++
 src/ModelOrderReduction.jl |   5 +-
 test/SOD.jl                | 139 +++++++++++++++++++++
 test/runtests.jl           |   3 +
 4 files changed, 388 insertions(+), 2 deletions(-)
 create mode 100644 src/DataReduction/SOD.jl
 create mode 100644 test/SOD.jl

diff --git a/src/DataReduction/SOD.jl b/src/DataReduction/SOD.jl
new file mode 100644
index 0000000..b99fb19
--- /dev/null
+++ b/src/DataReduction/SOD.jl
@@ -0,0 +1,243 @@
+"""
+    SOD <: AbstractDRProblem
+
+Smooth Orthogonal Decomposition (SOD) for extracting linear normal modes and natural
+frequencies from time-series data.
+
+SOD extends Proper Orthogonal Decomposition (POD) by considering temporal smoothness
+in addition to spatial variance. It uses the covariance matrices of displacement and
+velocity responses to form a generalized eigenvalue problem, which yields modes and
+frequencies without requiring knowledge of the system's mass matrix.
+
+# Fields
+- `snapshots`: Displacement snapshot data (matrix or vector of vectors), size (n_dof, n_samples)
+- `velocities`: Velocity snapshot data (same format as snapshots), or `:auto` to compute from snapshots
+- `dt`: Time step for computing velocities when `velocities = :auto`
+- `min_renergy`: Minimum relative energy threshold (0.0 to 1.0)
+- `min_nmodes`: Minimum number of modes to retain
+- `max_nmodes`: Maximum number of modes to retain
+- `nmodes`: Number of modes (computed after reduction)
+- `rbasis`: Reduced basis matrix (smooth orthogonal modes)
+- `renergy`: Relative energy captured by the modes
+- `spectrum`: Smooth orthogonal values (SOVs), which approximate ω² (squared frequencies)
+- `frequencies`: Estimated natural frequencies (√spectrum)
+
+# Constructor
+    SOD(snapshots; velocities=:auto, dt=1.0, min_renergy=1.0, min_nmodes=1, max_nmodes=size(snapshots,1))
+    SOD(snapshots, nmodes; velocities=:auto, dt=1.0)
+
+# References
+- Chelidze D, Zhou W. Smooth orthogonal decomposition-based vibration mode identification.
+  Journal of Sound and Vibration. 2006; 292(3-5): 461-473.
+  https://www.sciencedirect.com/science/article/abs/pii/S0022460X05005948
+"""
+mutable struct SOD <: AbstractDRProblem
+    # specified
+    snapshots::Any
+    velocities::Any
+    dt::Float64
+    min_renergy::Any
+    min_nmodes::Int
+    max_nmodes::Int
+    # computed
+    nmodes::Int
+    rbasis::Any
+    renergy::Any
+    spectrum::Any
+    frequencies::Any
+    # constructors
+    function SOD(
+            snaps;
+            velocities = :auto,
+            dt::Real = 1.0,
+            min_renergy = 1.0,
+            min_nmodes::Int = 1,
+            max_nmodes::Int = _get_ndof(snaps)
+        )
+        nmodes = min_nmodes
+        errorhandle_sod(snaps, velocities, dt, nmodes, min_renergy, min_nmodes, max_nmodes)
+        return new(
+            snaps, velocities, Float64(dt), min_renergy, min_nmodes, max_nmodes,
+            nmodes, missing, 1.0, missing, missing
+        )
+    end
+    function SOD(snaps, nmodes::Int; velocities = :auto, dt::Real = 1.0)
+        errorhandle_sod(snaps, velocities, dt, nmodes, 0.0, nmodes, nmodes)
+        return new(
+            snaps, velocities, Float64(dt), 0.0, nmodes, nmodes,
+            nmodes, missing, 1.0, missing, missing
+        )
+    end
+end
+
+# Helper to get number of degrees of freedom from snapshots
+function _get_ndof(snaps::AbstractMatrix)
+    return size(snaps, 1)
+end
+
+function _get_ndof(snaps::AbstractVector{<:AbstractVector})
+    return length(snaps[1])
+end
+
+# Error handling for SOD
+function errorhandle_sod(
+        data::AbstractMatrix, velocities, dt, nmodes::Int, max_energy,
+        min_nmodes, max_nmodes
+    )
+    @assert size(data, 1) > 1 "State vector is expected to be vector valued."
+    s = minimum(size(data))
+    @assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
+    @assert dt > 0 "Time step dt must be positive"
+    if velocities !== :auto
+        @assert size(velocities) == size(data) "Velocities must have same size as snapshots"
+    end
+    return nothing
+end
+
+function errorhandle_sod(
+        data::AbstractVector{T}, velocities, dt, nmodes::Int, max_energy, min_nmodes,
+        max_nmodes
+    ) where {T <: AbstractVector}
+    @assert size(data[1], 1) > 1 "State vector is expected to be vector valued."
+    s = min(size(data, 1), size(data[1], 1))
+    @assert 0 < nmodes <= s "Number of modes should be in {1,2,...,$s}."
+    @assert min_nmodes <= max_nmodes "Minimum number of modes must lie below maximum number of modes"
+    @assert 0.0 <= max_energy <= 1.0 "Maximum relative energy must be in [0,1]"
+    @assert dt > 0 "Time step dt must be positive"
+    if velocities !== :auto
+        @assert length(velocities) == length(data) "Velocities must have same length as snapshots"
+    end
+    return nothing
+end
+
+# Compute velocity from displacement using finite differences
+function _compute_velocity(X::AbstractMatrix, dt::Real)
+    # Central difference for interior points, forward/backward for boundaries
+    n_dof, n_samples = size(X)
+    V = similar(X)
+
+    if n_samples == 1
+        V .= 0
+        return V
+    elseif n_samples == 2
+        # Forward difference only
+        V[:, 1] = (X[:, 2] - X[:, 1]) / dt
+        V[:, 2] = V[:, 1]
+        return V
+    end
+
+    # Forward difference for first point
+    V[:, 1] = (X[:, 2] - X[:, 1]) / dt
+
+    # Central difference for interior points
+    for i in 2:(n_samples - 1)
+        V[:, i] = (X[:, i + 1] - X[:, i - 1]) / (2 * dt)
+    end
+
+    # Backward difference for last point
+    V[:, n_samples] = (X[:, n_samples] - X[:, n_samples - 1]) / dt
+
+    return V
+end
+
+function _compute_velocity(X::AbstractVector{<:AbstractVector}, dt::Real)
+    X_mat = matricize(X)
+    V_mat = _compute_velocity(X_mat, dt)
+    return V_mat
+end
+
+# Convert velocities to matrix form
+function _get_velocity_matrix(sod::SOD)
+    if sod.velocities === :auto
+        X = sod.snapshots isa AbstractMatrix ? sod.snapshots : matricize(sod.snapshots)
+        return _compute_velocity(X, sod.dt)
+    else
+        return sod.velocities isa AbstractMatrix ? sod.velocities : matricize(sod.velocities)
+    end
+end
+
+# Get displacement matrix
+function _get_displacement_matrix(sod::SOD)
+    return sod.snapshots isa AbstractMatrix ? sod.snapshots : matricize(sod.snapshots)
+end
+
+"""
+    reduce!(sod::SOD)
+
+Perform Smooth Orthogonal Decomposition on the snapshot data.
+
+This solves the generalized eigenvalue problem:
+    Σ_xx Ψ = λ Σ_vv Ψ
+
+where Σ_xx = X X^T / N (displacement covariance) and Σ_vv = V V^T / N (velocity covariance).
+
+The eigenvalues λ (smooth orthogonal values, SOVs) approximate ω² (squared natural frequencies).
+The eigenvectors Ψ are the smooth orthogonal modes (SOMs).
+"""
+function reduce!(sod::SOD)
+    X = _get_displacement_matrix(sod)
+    V = _get_velocity_matrix(sod)
+
+    n_dof, n_samples = size(X)
+
+    # Compute covariance matrices
+    # Σ_xx = X * X' / n_samples (displacement covariance)
+    # Σ_vv = V * V' / n_samples (velocity covariance)
+    Σ_xx = X * X' / n_samples
+    Σ_vv = V * V' / n_samples
+
+    # Add small regularization to Σ_vv if it's singular
+    # This handles cases where velocity variance is very small
+    eps_reg = 1.0e-10 * maximum(abs, Σ_vv)
+    if eps_reg < 1.0e-14
+        eps_reg = 1.0e-14
+    end
+    Σ_vv_reg = Σ_vv + eps_reg * I
+
+    # Solve generalized eigenvalue problem: Σ_xx Ψ = λ Σ_vv Ψ
+    # This is equivalent to finding eigenvalues of Σ_vv^(-1) Σ_xx
+    # For numerical stability, we use eigen on the symmetric pencil
+    eig_result = eigen(Symmetric(Σ_xx), Symmetric(Σ_vv_reg))
+
+    # Eigenvalues are the SOVs (smooth orthogonal values)
+    # They should be sorted in decreasing order (largest eigenvalue = highest energy)
+    sovs = real.(eig_result.values)
+    soms = real.(eig_result.vectors)
+
+    # Sort by eigenvalue magnitude (largest first)
+    perm = sortperm(sovs, rev = true)
+    sovs = sovs[perm]
+    soms = soms[:, perm]
+
+    # Determine truncation based on energy criterion
+    sod.nmodes, sod.renergy = determine_truncation(
+        sovs, sod.min_nmodes, sod.min_renergy, sod.max_nmodes
+    )
+
+    # Store results
+    sod.rbasis = soms[:, 1:(sod.nmodes)]
+    sod.spectrum = sovs
+
+    # Compute natural frequencies from SOVs
+    # SOVs approximate ω², so frequencies are √SOVs
+    # Only take sqrt of positive values
+    sod.frequencies = sqrt.(max.(sovs[1:(sod.nmodes)], 0.0))
+
+    return nothing
+end
+
+function Base.show(io::IO, sod::SOD)
+    print(io, "SOD \n")
+    print(io, "Reduction Order = ", sod.nmodes, "\n")
+    n_dof = _get_ndof(sod.snapshots)
+    n_samples = sod.snapshots isa AbstractMatrix ? size(sod.snapshots, 2) :
+        length(sod.snapshots)
+    print(io, "Snapshot size = (", n_dof, ",", n_samples, ")\n")
+    print(io, "Relative Energy = ", sod.renergy, "\n")
+    if !ismissing(sod.frequencies)
+        print(io, "Estimated frequencies = ", sod.frequencies, "\n")
+    end
+    return nothing
+end
diff --git a/src/ModelOrderReduction.jl b/src/ModelOrderReduction.jl
index 270084e..d06eae2 100644
--- a/src/ModelOrderReduction.jl
+++ b/src/ModelOrderReduction.jl
@@ -5,7 +5,7 @@ using DocStringExtensions: DocStringExtensions, FUNCTIONNAME, SIGNATURES, TYPEDS
 using ModelingToolkit: ModelingToolkit, @variables, Differential, Equation, Num, ODESystem,
     SymbolicUtils, Symbolics, arguments, build_function, complete,
     expand, substitute, tearing_substitution
-using LinearAlgebra: LinearAlgebra, /, \, mul!, qr, svd
+using LinearAlgebra: LinearAlgebra, /, \, I, Symmetric, eigen, mul!, qr, svd
 
 using Setfield: Setfield, @set!
 
@@ -15,8 +15,9 @@ include("Types.jl")
 include("ErrorHandle.jl")
 
 include("DataReduction/POD.jl")
+include("DataReduction/SOD.jl")
 export SVD, TSVD, RSVD
-export POD, reduce!
+export POD, SOD, reduce!
 
 include("deim.jl")
 export deim
diff --git a/test/SOD.jl b/test/SOD.jl
new file mode 100644
index 0000000..a4e2c54
--- /dev/null
+++ b/test/SOD.jl
@@ -0,0 +1,139 @@
+using Test, ModelOrderReduction
+using OrdinaryDiffEq
+
+# Simple harmonic oscillator system for testing SOD
+# This system has known natural frequencies that SOD should recover
+function harmonic_oscillator_prob()
+    # 2-DOF mass-spring system with known frequencies
+    # m1*x1'' + k1*x1 + k2*(x1-x2) = 0
+    # m2*x2'' + k2*(x2-x1) = 0
+    # With m1=m2=1, k1=1, k2=1, natural frequencies are ω₁ ≈ 0.618, ω₂ ≈ 1.618
+
+    function harmonic!(du, u, p, t)
+        # u = [x1, v1, x2, v2]
+        k1, k2, m1, m2 = p
+        x1, v1, x2, v2 = u
+        du[1] = v1
+        du[2] = (-k1 * x1 - k2 * (x1 - x2)) / m1
+        du[3] = v2
+        return du[4] = (-k2 * (x2 - x1)) / m2
+    end
+
+    u0 = [1.0, 0.0, 0.5, 0.0]  # Initial displacement, zero velocity
+    p = [1.0, 1.0, 1.0, 1.0]   # k1, k2, m1, m2
+    tspan = (0.0, 50.0)
+    prob = ODEProblem(harmonic!, u0, tspan, p)
+    sol = solve(prob, Tsit5(), saveat = 0.01)
+    return sol
+end
+
+# Use Lorenz system for comparison with POD (same test as POD uses)
+function lorenz_prob()
+    function lorenz!(du, u, p, t)
+        du[1] = p[1] * (u[2] - u[1])
+        du[2] = u[1] * (p[2] - u[3]) - u[2]
+        return du[3] = u[1] * u[2] - p[3] * u[3]
+    end
+
+    u0 = [1, 0, 0]
+    p = [10, 28, 8 / 3]
+    tspan = (0, 100)
+    prob = ODEProblem(lorenz!, u0, tspan, p)
+    sol = solve(prob, Tsit5(), saveat = 0.1)
+    return sol
+end
+
+@testset "SOD - Basic Functionality" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Test with matrix input
+    order = 2
+    matrix_reducer = SOD(solution, order; dt = dt)
+    reduce!(matrix_reducer)
+
+    @test size(matrix_reducer.rbasis, 2) == matrix_reducer.nmodes
+    @test size(matrix_reducer.rbasis, 1) == size(solution, 1)
+    @test !ismissing(matrix_reducer.spectrum)
+    @test !ismissing(matrix_reducer.frequencies)
+    @test length(matrix_reducer.frequencies) == order
+
+    # Test with vector of vectors input
+    snapshot_reducer = SOD(sol.u, order; dt = dt)
+    reduce!(snapshot_reducer)
+
+    @test size(snapshot_reducer.rbasis, 2) == snapshot_reducer.nmodes
+    @test !ismissing(snapshot_reducer.frequencies)
+end
+
+@testset "SOD - Energy Truncation" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Test with energy-based truncation
+    reducer = SOD(solution; dt = dt, min_nmodes = 1, max_nmodes = 3, min_renergy = 0.1)
+    reduce!(reducer)
+    @test reducer.renergy >= 0.1
+    @test 1 <= reducer.nmodes <= 3
+end
+
+@testset "SOD - Provided Velocities" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    # Manually compute velocities (central difference)
+    n_dof, n_samples = size(solution)
+    velocities = similar(solution)
+    velocities[:, 1] = (solution[:, 2] - solution[:, 1]) / dt
+    for i in 2:(n_samples - 1)
+        velocities[:, i] = (solution[:, i + 1] - solution[:, i - 1]) / (2 * dt)
+    end
+    velocities[:, n_samples] = (solution[:, n_samples] - solution[:, n_samples - 1]) / dt
+
+    # Test with provided velocities
+    order = 2
+    reducer = SOD(solution, order; velocities = velocities, dt = dt)
+    reduce!(reducer)
+
+    @test size(reducer.rbasis, 2) == reducer.nmodes
+    @test !ismissing(reducer.frequencies)
+end
+
+@testset "SOD - Harmonic Oscillator Frequency Recovery" begin
+    sol = harmonic_oscillator_prob()
+
+    # Extract position data (x1 and x2, not velocities)
+    positions = vcat(sol[1, :]', sol[3, :]')  # Shape: (2, n_samples)
+    dt = 0.01
+
+    # SOD should find the modes of this oscillating system
+    reducer = SOD(positions, 2; dt = dt)
+    reduce!(reducer)
+
+    @test reducer.nmodes == 2
+    @test !ismissing(reducer.frequencies)
+
+    # Check that we get reasonable positive frequencies
+    # (The exact values depend on the sampling and system dynamics)
+    @test all(reducer.frequencies .>= 0)
+end
+
+@testset "SOD - Show Method" begin
+    sol = lorenz_prob()
+    solution = Array(sol)
+    dt = 0.1
+
+    reducer = SOD(solution, 2; dt = dt)
+    reduce!(reducer)
+
+    # Test that show doesn't error
+    io = IOBuffer()
+    show(io, reducer)
+    str = String(take!(io))
+    @test occursin("SOD", str)
+    @test occursin("Reduction Order", str)
+    @test occursin("frequencies", str)
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 6d1c1af..0e6e978 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,9 @@ end
 @safetestset "POD" begin
     include("DataReduction.jl")
 end
+@safetestset "SOD" begin
+    include("SOD.jl")
+end
 @safetestset "utils" begin
     include("utils.jl")
 end