""" Chaos analysis: Lyapunov spectrum, multifractal formalism, RQA. Implements comprehensive chaos metrics including Lyapunov exponents, Kolmogorov-Sinai entropy, correlation dimension, multifractal spectrum, and recurrence quantification analysis. """ import numpy as np from scipy.stats import entropy from scipy.special import factorial from typing import Dict, List, Optional class ChaosAnalyzer: """ Comprehensive chaos prediction and analysis system. """ def __init__(self, r: float = 3.9, x0: float = 0.4): self.r = r self.state = x0 self.history: List[float] = [] self.lyapunov_spectrum: List[float] = [] self.kolmogorov_sinai_entropy: Optional[float] = None self.correlation_dim: Optional[float] = None self.generalized_dimensions: Dict[float, float] = {} self.multifractal_spectrum: Dict[str, list] = {} self.recurrence_quantification: Dict[str, float] = {} def run_logistic_map(self, steps: int): """Run logistic map evolution.""" for _ in range(steps): self.state = self.r * self.state * (1 - self.state) self.history.append(self.state) def complete_analysis(self, steps: int): """Perform comprehensive chaos analysis.""" self.run_logistic_map(steps) self.calculate_lyapunov_spectrum() self.calculate_kolmogorov_sinai_entropy() self.calculate_correlation_dimension() self.calculate_information_metrics() self.calculate_multifractal_spectrum() self.recurrence_quantification_analysis() def calculate_lyapunov_spectrum(self) -> Optional[float]: """Lyapunov exponent via Jacobian method.""" if len(self.history) < 100: return None lyapunov_sum = 0.0 count = 0 for x in self.history[50:]: if 0 < x < 1: derivative = abs(self.r * (1 - 2 * x)) if derivative > 0: lyapunov_sum += np.log(derivative) count += 1 if count == 0: return 0.0 exp = lyapunov_sum / count self.lyapunov_spectrum = [exp] return exp def calculate_kolmogorov_sinai_entropy(self) -> float: """KS entropy: h_KS = sum(lambda_i for lambda_i > 0).""" positive = [l for l in self.lyapunov_spectrum if l > 0] self.kolmogorov_sinai_entropy = sum(positive) return self.kolmogorov_sinai_entropy def calculate_correlation_dimension(self) -> Optional[float]: """Correlation dimension via Grassberger-Procaccia algorithm.""" if len(self.history) < 100: return None data = np.array(self.history[-1000:]) m, tau = 2, 1 N = len(data) - (m - 1) * tau embedded = np.zeros((N, m)) for i in range(m): embedded[:, i] = data[i * tau : i * tau + N] radii = np.logspace(-3, 0, 20) C_r = [] for r in radii: count = 0 for i in range(N): dists = np.linalg.norm(embedded[i] - embedded[i + 1 :], axis=1) count += np.sum(dists < r) C_r.append(2 * count / (N * (N - 1)) if N > 1 else 0) C_r = np.array(C_r) valid = C_r > 0 if np.sum(valid) > 2: slope, _ = np.polyfit(np.log(radii[valid]), np.log(C_r[valid]), 1) self.correlation_dim = slope else: self.correlation_dim = 0.0 return self.correlation_dim def calculate_information_metrics(self) -> Dict[str, float]: """Mutual information and permutation entropy.""" if len(self.history) < 100: return {} data = np.array(self.history) results = {} # Mutual information at various delays for delay in [1, 5, 10]: if len(data) > delay: mi = self._mutual_information(data[:-delay], data[delay:]) results[f'mutual_info_lag_{delay}'] = mi # Permutation entropy results['permutation_entropy'] = self._permutation_entropy(data) return results def _mutual_information(self, x: np.ndarray, y: np.ndarray) -> float: """Mutual information between two series.""" bins = 20 hist_2d, _, _ = np.histogram2d(x, y, bins=bins) p_xy = hist_2d / np.sum(hist_2d) p_x = np.sum(p_xy, axis=1) p_y = np.sum(p_xy, axis=0) mi = 0.0 for i in range(bins): for j in range(bins): if p_xy[i, j] > 0 and p_x[i] > 0 and p_y[j] > 0: mi += p_xy[i, j] * np.log(p_xy[i, j] / (p_x[i] * p_y[j])) return mi def _permutation_entropy(self, data: np.ndarray, order: int = 3) -> float: """Permutation entropy of ordinal patterns.""" if len(data) < order: return 0.0 patterns = [] for i in range(len(data) - order + 1): pattern = tuple(np.argsort(data[i : i + order])) patterns.append(pattern) _, counts = np.unique(patterns, axis=0, return_counts=True) probs = counts / len(patterns) pe = entropy(probs, base=2) max_ent = np.log2(factorial(order, exact=True)) return pe / max_ent if max_ent > 0 else 0.0 def calculate_multifractal_spectrum(self): """Multifractal spectrum via generalized dimensions D_q.""" if len(self.history) < 500: return data = np.array(self.history) q_values = np.linspace(-5, 5, 21) boxes = 20 hist, _ = np.histogram(data, bins=boxes) probs = hist / np.sum(hist) probs = probs[probs > 0] for q in q_values: if q == 1: D_q = entropy(probs, base=2) / np.log(boxes) else: Z_q = np.sum(probs ** q) D_q = np.log(Z_q) / ((1 - q) * np.log(boxes)) if Z_q > 0 else 0 self.generalized_dimensions[q] = D_q # f(alpha) spectrum via Legendre transform qs = np.array(list(self.generalized_dimensions.keys())) Dqs = np.array(list(self.generalized_dimensions.values())) alphas, f_alphas = [], [] for i in range(1, len(qs) - 1): alpha = (Dqs[i + 1] - Dqs[i - 1]) / (qs[i + 1] - qs[i - 1]) f_alpha = qs[i] * alpha - qs[i] * Dqs[i] alphas.append(alpha) f_alphas.append(f_alpha) self.multifractal_spectrum = { 'alpha': alphas, 'f_alpha': f_alphas, 'width': max(alphas) - min(alphas) if alphas else 0, } def recurrence_quantification_analysis(self): """RQA: recurrence rate, determinism, laminarity, entropy.""" if len(self.history) < 100: return data = np.array(self.history[-500:]) N = len(data) threshold = 0.1 * np.std(data) # Vectorized recurrence matrix dists = np.abs(data[:, None] - data[None, :]) R = (dists < threshold).astype(float) RR = np.sum(R) / (N * N) # Determinism (diagonal line ratio) diag_sums = [] for k in range(1, N): diag = np.diag(R, k) s = np.sum(diag) if s > 2: diag_sums.append(s) total_R = np.sum(R) DET = sum(diag_sums) / total_R if total_R > 0 else 0 # Entropy of diagonal line lengths if diag_sums: diag_probs = np.array(diag_sums) / sum(diag_sums) ENTR = entropy(diag_probs, base=2) else: ENTR = 0.0 self.recurrence_quantification = { 'recurrence_rate': RR, 'determinism': DET, 'entropy': ENTR, 'is_chaotic': DET < 0.99 and RR > 0.01, }