"""Model-data compatibility checking for rheological models.
This module provides functions to assess whether a given model is appropriate
for a dataset based on the underlying physics and data characteristics.
The compatibility checker helps users understand when model failures are expected
due to physics mismatch rather than optimization issues.
"""
from __future__ import annotations
from enum import Enum
from typing import TYPE_CHECKING, Any
import numpy as np
from rheojax.logging import get_logger
if TYPE_CHECKING:
from rheojax.core.base import BaseModel
logger = get_logger(__name__)
[docs]
class DecayType(Enum):
"""Types of relaxation decay behavior."""
EXPONENTIAL = "exponential" # Simple Maxwell-like exp(-t/tau)
POWER_LAW = "power_law" # Power-law t^(-alpha)
STRETCHED = "stretched" # Stretched exponential exp(-(t/tau)^beta)
MITTAG_LEFFLER = "mittag_leffler" # Mittag-Leffler E_alpha(-(t/tau)^alpha)
MULTI_MODE = "multi_mode" # Multiple relaxation modes
UNKNOWN = "unknown" # Cannot determine
[docs]
class MaterialType(Enum):
"""Types of material behavior."""
SOLID = "solid" # Solid-like (finite equilibrium modulus)
LIQUID = "liquid" # Liquid-like (zero equilibrium modulus, flows)
GEL = "gel" # Gel-like (power-law relaxation)
VISCOELASTIC_SOLID = "viscoelastic_solid" # Viscoelastic solid
VISCOELASTIC_LIQUID = "viscoelastic_liquid" # Viscoelastic liquid
UNKNOWN = "unknown"
[docs]
def detect_decay_type(t: np.ndarray, G_t: np.ndarray) -> DecayType:
"""Detect the type of relaxation decay from time-domain data.
Analyzes the decay pattern to determine if it follows exponential,
power-law, stretched exponential, or Mittag-Leffler behavior.
Parameters
----------
t : np.ndarray
Time array (s)
G_t : np.ndarray
Relaxation modulus array (Pa)
Returns
-------
DecayType
Detected decay type
"""
logger.debug(
"Detecting decay type",
n_points=len(t),
t_range=(float(t.min()), float(t.max())) if len(t) > 0 else (0, 0),
)
from scipy.stats import linregress
if len(t) < 10 or len(G_t) < 10:
logger.debug("Insufficient data points for decay detection", n_points=len(t))
return DecayType.UNKNOWN
# Remove any invalid values
valid = np.isfinite(t) & np.isfinite(G_t) & (t > 0) & (G_t > 0)
if np.sum(valid) < 10:
n_non_positive = int(np.sum(np.asarray(G_t) <= 0))
if n_non_positive > 0:
logger.info(
"Cannot classify decay type: G(t) contains non-positive values",
n_non_positive=n_non_positive,
data_range=(float(np.min(G_t)), float(np.max(G_t))),
)
else:
logger.debug(
"Insufficient valid data points after filtering",
n_valid=int(np.sum(valid)),
)
return DecayType.UNKNOWN
t = t[valid]
G_t = G_t[valid]
# Ensure data is sorted by time
sort_idx = np.argsort(t)
t = t[sort_idx]
G_t = G_t[sort_idx]
# R8-COMPAT-001: narrow time windows give unreliable decay classification.
# Threshold: 0.5 decades (~3x time range) catches near-constant data
# while preserving classification for short-duration industrial QC tests.
t_span_decades = np.log10(t.max() / max(t.min(), 1e-30))
if t_span_decades < 0.5:
logger.info("Time span < 0.5 decades; decay type classification unreliable")
return DecayType.UNKNOWN
# 1. Check for exponential decay: log(G) vs t should be linear
# Normalize t to [0,1] so R² is comparable with the log-log power-law R²
slope_exp = 0.0
try:
with np.errstate(divide="ignore", invalid="ignore"):
log_G = np.log(np.maximum(G_t, 1e-30))
t_range = t.max() - t.min()
t_norm = (t - t.min()) / (t_range + 1e-30) if t_range > 0 else t
slope_exp, intercept_exp, r_exp, _, _ = linregress(t_norm, log_G)
r_exp_sq = float(r_exp**2) if np.isfinite(r_exp) else 0.0
logger.debug(
"Exponential fit",
r_squared=float(r_exp_sq),
slope=float(slope_exp),
)
except ValueError as e:
r_exp_sq = 0.0
logger.debug("Exponential fit failed", error=str(e))
# 2. Check for power-law decay: log(G) vs log(t) should be linear
slope_pow = 0.0
try:
with np.errstate(divide="ignore", invalid="ignore"):
log_t = np.log(np.maximum(t, 1e-30))
log_G = np.log(np.maximum(G_t, 1e-30))
slope_pow, intercept_pow, r_pow, _, _ = linregress(log_t, log_G)
r_pow_sq = float(r_pow**2) if np.isfinite(r_pow) else 0.0
logger.debug(
"Power-law fit",
r_squared=float(r_pow_sq),
slope=float(slope_pow),
)
except ValueError as e:
r_pow_sq = 0.0
logger.debug("Power-law fit failed", error=str(e))
# 3. Check for stretched exponential: log(-log(G/G0)) vs log(t)
# COMPAT-001: Initialize slope_stretch before the try block so the variable
# is always defined at the decision logic below. Python's short-circuit
# evaluation means slope_stretch is never used when r_stretch_sq == 0.0,
# but initializing explicitly eliminates the UnboundLocalError risk if
# threshold values change in the future.
slope_stretch = 0.0
try:
G0 = G_t[0]
G_norm = G_t / G0
G_norm = np.clip(G_norm, 1e-10, 1.0) # Avoid log(0)
# Suppress expected warnings for invalid log operations (handled by isfinite check)
with np.errstate(divide="ignore", invalid="ignore"):
log_log = np.log(-np.log(G_norm))
log_t = np.log(t)
valid_stretch = np.isfinite(log_log)
if np.sum(valid_stretch) >= 5:
slope_stretch, _, r_stretch, _, _ = linregress(
log_t[valid_stretch], log_log[valid_stretch]
)
r_stretch_sq = float(r_stretch**2) if np.isfinite(r_stretch) else 0.0
logger.debug(
"Stretched exponential fit",
r_squared=float(r_stretch_sq),
n_valid_points=int(np.sum(valid_stretch)),
)
else:
r_stretch_sq = 0.0
logger.debug("Stretched exponential fit: insufficient valid points")
except (ValueError, RuntimeWarning) as e:
r_stretch_sq = 0.0
logger.debug("Stretched exponential fit failed", error=str(e))
# Decision logic
threshold_high = 0.90 # High confidence
threshold_med = 0.75 # Medium confidence
# Exponential decay (Maxwell-like)
if r_exp_sq > threshold_high:
# Check if slope is negative (decay)
if slope_exp < 0:
logger.info(
"Detected exponential decay",
r_squared=float(r_exp_sq),
slope=float(slope_exp),
)
return DecayType.EXPONENTIAL
# Power-law decay (gel-like)
if r_pow_sq > threshold_high:
# Check if slope is negative (decay)
if slope_pow < 0:
logger.info(
"Detected power-law decay",
r_squared=float(r_pow_sq),
slope=float(slope_pow),
)
return DecayType.POWER_LAW
# Stretched exponential: R² must be high AND the slope (β exponent) must be ≤ 1.
# R10-COMPAT-001: without the slope guard, compressed exponentials (β > 1) are
# misclassified as STRETCHED because log(-log(G/G0)) vs log(t) is also linear for
# β > 1 — the only distinction is the slope magnitude.
if r_stretch_sq > threshold_high and 0 < slope_stretch <= 1.0:
logger.info(
"Detected stretched exponential decay",
r_squared=float(r_stretch_sq),
beta_exponent=float(slope_stretch),
)
return DecayType.STRETCHED
# Multi-mode if nothing fits well but data shows decay
if G_t[-1] < G_t[0] * 0.9: # At least 10% decay
# Check if it's a combination
if r_exp_sq > threshold_med or r_pow_sq > threshold_med:
logger.info(
"Detected multi-mode decay",
r_exp_sq=float(r_exp_sq),
r_pow_sq=float(r_pow_sq),
)
return DecayType.MULTI_MODE
else:
# Complex decay pattern, likely Mittag-Leffler
logger.info(
"Detected Mittag-Leffler-like decay",
decay_ratio=float(G_t[-1] / G_t[0]),
)
return DecayType.MITTAG_LEFFLER
logger.debug("Could not determine decay type")
return DecayType.UNKNOWN
[docs]
def detect_material_type(
t: np.ndarray | None = None,
G_t: np.ndarray | None = None,
omega: np.ndarray | None = None,
G_star: np.ndarray | None = None,
) -> MaterialType:
"""Detect the material type from relaxation or oscillation data.
Parameters
----------
t : np.ndarray, optional
Time array for relaxation data (s)
G_t : np.ndarray, optional
Relaxation modulus (Pa)
omega : np.ndarray, optional
Frequency array for oscillation data (rad/s)
G_star : np.ndarray, optional
Complex modulus array with shape (N, 2) where [:, 0] is G' and [:, 1] is G"
Returns
-------
MaterialType
Detected material type
"""
# Use relaxation data if available
if t is not None and G_t is not None and len(t) >= 10:
return _detect_from_relaxation(t, G_t)
# Use oscillation data
if omega is not None and G_star is not None and len(omega) >= 10:
return _detect_from_oscillation(omega, G_star)
return MaterialType.UNKNOWN
def _detect_from_relaxation(t: np.ndarray, G_t: np.ndarray) -> MaterialType:
"""Detect material type from relaxation modulus."""
# Check for equilibrium modulus (finite G at long times)
if len(G_t) < 10:
return MaterialType.UNKNOWN
# Check decay type first - power-law is special
decay_type = detect_decay_type(t, G_t)
# Power-law relaxation is characteristic of gels, regardless of decay magnitude
if decay_type == DecayType.POWER_LAW:
# Power-law materials are gels or viscoelastic solids
# Distinguish based on equilibrium modulus
G_final = np.mean(G_t[-5:]) # Average last 5 points
G_initial = G_t[0]
decay_ratio = G_final / G_initial if G_initial != 0 else 0.0
if decay_ratio > 0.3: # Significant final modulus
return MaterialType.VISCOELASTIC_SOLID
else:
return MaterialType.GEL
# For non-power-law materials, use decay ratio
G_final = np.mean(G_t[-5:]) # Average last 5 points
G_initial = np.mean(G_t[:5]) if len(G_t) >= 5 else G_t[0]
# R8-COMPAT-002: guard against zero G_initial
decay_ratio = G_final / G_initial if G_initial != 0 else 0.0
if decay_ratio > 0.5:
# Finite equilibrium modulus → solid-like
return MaterialType.VISCOELASTIC_SOLID
elif decay_ratio < 0.1:
# Strong decay → liquid-like
return MaterialType.VISCOELASTIC_LIQUID
else:
# Intermediate behavior
return MaterialType.UNKNOWN
def _detect_from_oscillation(omega: np.ndarray, G_star: np.ndarray) -> MaterialType:
"""Detect material type from complex modulus."""
# SUP-002: Handle 1D and complex inputs gracefully
if np.iscomplexobj(G_star):
# Convert complex to [real, imag] column stack
G_star = np.column_stack([np.real(G_star), np.imag(G_star)])
if G_star.ndim != 2 or G_star.shape[1] != 2:
return MaterialType.UNKNOWN
G_prime = G_star[:, 0] # Storage modulus
G_double_prime = G_star[:, 1] # Loss modulus
# Low-frequency behavior
low_freq_idx = np.argmin(omega)
# Check if G' > G" at low frequency (solid-like)
# or G" > G' (liquid-like)
if G_prime[low_freq_idx] > G_double_prime[low_freq_idx]:
return MaterialType.SOLID
elif G_double_prime[low_freq_idx] > G_prime[low_freq_idx]:
return MaterialType.LIQUID
else:
return MaterialType.UNKNOWN
[docs]
def check_model_compatibility(
model: BaseModel,
t: np.ndarray | None = None,
G_t: np.ndarray | None = None,
omega: np.ndarray | None = None,
G_star: np.ndarray | None = None,
test_mode: str | None = None,
) -> dict[str, Any]:
"""Check if a model is compatible with the given data.
This function analyzes the data characteristics and compares them with
the model's underlying physics to assess compatibility.
Parameters
----------
model : BaseModel
The rheological model to check
t : np.ndarray, optional
Time array for relaxation data (s)
G_t : np.ndarray, optional
Relaxation modulus (Pa)
omega : np.ndarray, optional
Frequency array for oscillation data (rad/s)
G_star : np.ndarray, optional
Complex modulus array with shape (N, 2)
test_mode : str, optional
Test mode ('relaxation', 'creep', 'oscillation')
Returns
-------
dict
Dictionary with compatibility information:
- 'compatible': bool, whether model is likely compatible
- 'confidence': float, confidence level (0-1)
- 'decay_type': DecayType, detected decay pattern
- 'material_type': MaterialType, detected material behavior
- 'warnings': list[str], compatibility warnings
- 'recommendations': list[str], suggested alternative models
"""
warnings = []
recommendations = []
# Detect data characteristics
decay_type = DecayType.UNKNOWN
material_type = MaterialType.UNKNOWN
if test_mode == "relaxation" and t is not None and G_t is not None:
decay_type = detect_decay_type(t, G_t)
material_type = detect_material_type(t=t, G_t=G_t)
elif test_mode == "oscillation" and omega is not None and G_star is not None:
material_type = detect_material_type(omega=omega, G_star=G_star)
elif test_mode is not None and test_mode not in ("relaxation", "oscillation"):
# COMPAT-002: Compatibility checking only supports relaxation and
# oscillation data. For other test_modes (flow_curve, startup, creep,
# laos), we return an uninformative result. Log a warning so callers
# are not silently misled into thinking a full compatibility check ran.
logger.warning(
"check_model_compatibility: test_mode=%r is not supported for "
"compatibility analysis (only 'relaxation' and 'oscillation' are "
"implemented). Returning default compatible=True with low confidence.",
test_mode,
)
# R7-COMPAT-001: Return early for unsupported test modes.
# Previously, confidence=0.3 set here was overwritten by the default
# confidence=0.5 below, silently discarding the low-confidence signal.
return {
"compatible": True,
"confidence": 0.3,
"decay_type": decay_type,
"material_type": material_type,
"warnings": [
f"Compatibility checking not implemented for test_mode={test_mode!r}. "
"Only 'relaxation' and 'oscillation' are supported."
],
"recommendations": [],
}
# Get model name
model_name = model.__class__.__name__
# Define model requirements
compatible = True
confidence = 0.5 # Default medium confidence
# Fractional Zener Solid-Solid (FZSS)
if "FractionalZenerSolidSolid" in model_name or model_name == "FZSS":
# FZSS expects Mittag-Leffler decay with finite equilibrium modulus
if decay_type == DecayType.EXPONENTIAL:
compatible = False
confidence = 0.9
warnings.append(
"FZSS model expects Mittag-Leffler (power-law) relaxation, "
"but data shows exponential decay."
)
recommendations.append("Maxwell")
recommendations.append("Zener")
elif decay_type == DecayType.POWER_LAW:
compatible = True
confidence = 0.8
elif material_type == MaterialType.VISCOELASTIC_LIQUID:
compatible = False
confidence = 0.7
warnings.append(
"FZSS model is designed for solid-like materials, "
"but data shows liquid-like behavior."
)
recommendations.append("FractionalMaxwellLiquid")
# Fractional Maxwell Liquid (FML)
elif "FractionalMaxwellLiquid" in model_name or model_name == "FML":
# FML expects liquid-like behavior (no equilibrium modulus)
if material_type == MaterialType.SOLID:
compatible = False
confidence = 0.8
warnings.append(
"FractionalMaxwellLiquid expects liquid-like behavior, "
"but data shows solid-like characteristics."
)
recommendations.append("FractionalZenerSolidSolid")
recommendations.append("FractionalKelvinVoigt")
# Fractional Maxwell Gel (FMG)
elif "FractionalMaxwellGel" in model_name or model_name == "FMG":
# FMG expects power-law relaxation (gel-like)
if decay_type == DecayType.EXPONENTIAL:
compatible = False
confidence = 0.85
warnings.append(
"FractionalMaxwellGel expects power-law relaxation, "
"but data shows exponential decay."
)
recommendations.append("Maxwell")
# Classical Maxwell model
elif model_name == "Maxwell":
# Maxwell expects exponential decay
if decay_type == DecayType.POWER_LAW:
compatible = False
confidence = 0.85
warnings.append(
"Maxwell model expects exponential decay, "
"but data shows power-law behavior."
)
recommendations.append("FractionalMaxwellGel")
recommendations.append("FractionalZenerSolidSolid")
# Classical Zener model
elif model_name == "Zener":
# Zener expects exponential decay with equilibrium modulus
if decay_type == DecayType.POWER_LAW:
compatible = False
confidence = 0.8
warnings.append(
"Zener model expects exponential decay, "
"but data shows power-law behavior."
)
recommendations.append("FractionalZenerSolidSolid")
# Fractional Kelvin-Voigt
elif "FractionalKelvinVoigt" in model_name and "Zener" not in model_name:
# FKV expects solid-like behavior
if material_type == MaterialType.VISCOELASTIC_LIQUID:
compatible = False
confidence = 0.75
warnings.append(
"FractionalKelvinVoigt is designed for solid-like materials, "
"but data shows liquid-like behavior."
)
recommendations.append("FractionalMaxwellLiquid")
return {
"compatible": compatible,
"confidence": confidence,
"decay_type": decay_type,
"material_type": material_type,
"warnings": warnings,
"recommendations": recommendations,
}
__all__ = [
"DecayType",
"MaterialType",
"detect_decay_type",
"detect_material_type",
"check_model_compatibility",
"format_compatibility_message",
]