#!/usr/bin/env python3 """pyGDIS 1.0.0 feature engineering and classification tutorial. Each training sample is an independent *parameter sweep*, not a single trajectory. Two synthetic populations have different hidden dynamics, but share the same observable parameter grid. Results are illustrative only. Run: python pygdis_ml_classification.py Outputs: ml_output/{features.csv,metrics.txt,*.png,*.pdf} """ from pathlib import Path import argparse import shutil import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from sklearn.model_selection import train_test_split, StratifiedKFold, cross_val_score from sklearn.pipeline import make_pipeline from sklearn.preprocessing import StandardScaler from sklearn.linear_model import LogisticRegression from sklearn.metrics import (balanced_accuracy_score, classification_report, ConfusionMatrixDisplay, RocCurveDisplay, roc_auc_score) from gdis import GDIS SEED = 2026 OUT = Path(__file__).resolve().parent / 'ml_output' OUT.mkdir(exist_ok=True) R = np.linspace(3.35, 3.79, 12) N_FAMILIES_PER_CLASS = 24 # 48 independent samples; increase for experimentation N_TRANSIENT, N_STEPS = 230, 220 def simulate_sweep(shift, rng): """One family: same measured r grid, different hidden r-offset. Tiny independent perturbations give distinct replicate families without encoding class in an explicit column. All r-effective remain below 4. """ trajectories = [] for r in R: effective_r = float(r + shift + rng.normal(0, 0.002)) x = float(rng.uniform(0.2, 0.8)) for _ in range(N_TRANSIENT): x = effective_r * x * (1 - x) values = np.empty(N_STEPS) for t in range(N_STEPS): x = effective_r * x * (1 - x) values[t] = x trajectories.append(values[:, None]) return trajectories def features_from_sweep(trajectories): """Fit pyGDIS to ONE independent sweep to prevent cross-family leakage.""" result = GDIS().fit_transform(trajectories=trajectories, parameters=R) score = np.asarray(result.gdis, dtype=float) sust = np.asarray(result.sustained_instability, dtype=float) energy = np.asarray(result.components['transition_energy'], dtype=float) all_states = np.concatenate([t[:, 0] for t in trajectories]) return { # A deliberately simple baseline (two aggregate time-series statistics) 'raw_mean': float(np.mean(all_states)), 'raw_std': float(np.std(all_states)), # Features of the parameter-ordered pyGDIS profile 'gdis_mean': float(np.mean(score)), 'gdis_std': float(np.std(score)), 'gdis_max': float(np.max(score)), 'gdis_range': float(np.ptp(score)), 'gdis_auc': float(np.trapezoid(score, R) / (R[-1] - R[0])) if hasattr(np, 'trapezoid') else float(np.trapz(score, R) / (R[-1] - R[0])), 'sustained_mean': float(np.mean(sust)), 'transition_energy_max': float(np.max(energy)), 'transition_energy_mean': float(np.mean(energy)), } def build_dataset(): rng = np.random.default_rng(SEED) records = [] for label, shift in ((0, -0.045), (1, +0.045)): for replicate in range(N_FAMILIES_PER_CLASS): sweep = simulate_sweep(shift, rng) row = features_from_sweep(sweep) row.update({'sample_id': f'class{label}_{replicate:03d}', 'label': label}) records.append(row) print(f"Processed {len(records):02d}/{2*N_FAMILIES_PER_CLASS} sweeps", end='\r', flush=True) print() df = pd.DataFrame(records) if not np.isfinite(df.drop(columns=['sample_id']).to_numpy()).all(): raise ValueError('Non-finite features. Inspect generated trajectories.') df.to_csv(OUT / 'features.csv', index=False) return df def evaluate(df): y = df['label'] train_idx, test_idx = train_test_split( np.arange(len(df)), test_size=0.30, random_state=SEED, stratify=y) baseline_cols = ['raw_mean', 'raw_std'] gdis_cols = [c for c in df.columns if c.startswith(('gdis_', 'sustained_', 'transition_'))] models = {} for name, columns in [('Raw baseline', baseline_cols), ('pyGDIS only', gdis_cols), ('Combined', baseline_cols + gdis_cols)]: model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000)) X_train = df.iloc[train_idx][columns] X_test = df.iloc[test_idx][columns] model.fit(X_train, y.iloc[train_idx]) prediction = model.predict(X_test) probability = model.predict_proba(X_test)[:, 1] models[name] = {'model': model, 'columns': columns, 'prediction': prediction, 'probability': probability, 'balanced_accuracy': balanced_accuracy_score(y.iloc[test_idx], prediction), 'auc': roc_auc_score(y.iloc[test_idx], probability)} return models, train_idx, test_idx def plot_figures(df, models, test_idx): plt.rcParams.update({'font.size': 11, 'axes.spines.top': False, 'axes.spines.right': False, 'savefig.facecolor': 'white'}) fig, axes = plt.subplots(1, 2, figsize=(11, 4.6), layout='constrained') for label, col, name in [(0, '#2563eb', 'Class 0'), (1, '#f97316', 'Class 1')]: d = df[df.label == label] axes[0].scatter(d.gdis_mean, d.transition_energy_max, label=name, c=col, s=58, alpha=0.82, edgecolors='white') axes[1].scatter(d.raw_mean, d.raw_std, label=name, c=col, s=58, alpha=0.82, edgecolors='white') axes[0].set(xlabel='Mean GDIS', ylabel='Maximum transition energy', title='Dynamical features') axes[1].set(xlabel='Mean state value', ylabel='State standard deviation', title='Raw statistical baseline') for ax in axes: ax.legend(frameon=False) ax.grid(alpha=.17) fig.savefig(OUT / '05_ml_feature_scatter.png', dpi=230) fig.savefig(OUT / '05_ml_feature_scatter.pdf') plt.close(fig) y_test = df.iloc[test_idx].label fig, axes = plt.subplots(1, 2, figsize=(11, 4.5), layout='constrained') for name, color in [('Raw baseline', '#64748b'), ('pyGDIS only', '#2563eb'), ('Combined', '#0d9488')]: RocCurveDisplay.from_predictions(y_test, models[name]['probability'], name=name, ax=axes[0]) axes[0].plot([0, 1], [0, 1], '--', color='#94a3b8') axes[0].set_title('Held-out ROC (one split)') ConfusionMatrixDisplay.from_predictions( y_test, models['Combined']['prediction'], ax=axes[1], cmap='Blues', colorbar=False, display_labels=['Class 0', 'Class 1']) axes[1].set_title('Combined model: held-out confusion matrix') fig.savefig(OUT / '06_ml_evaluation.png', dpi=230) fig.savefig(OUT / '06_ml_evaluation.pdf') plt.close(fig) names = list(models) fig, ax = plt.subplots(figsize=(9, 4.6), layout='constrained') bars = ax.bar(names, [models[n]['balanced_accuracy'] for n in names], color=['#94a3b8', '#2563eb', '#0d9488'], width=.58) ax.bar_label(bars, fmt='%.2f', padding=5) ax.set(ylabel='Held-out balanced accuracy', ylim=(0, 1.10), title='Baseline vs pyGDIS features (illustrative)') ax.grid(axis='y', alpha=.18) fig.savefig(OUT / '07_ml_model_comparison.png', dpi=230) fig.savefig(OUT / '07_ml_model_comparison.pdf') plt.close(fig) def main(): parser = argparse.ArgumentParser(description="Generate pyGDIS classification figures") parser.add_argument("--publish-static-dir", type=Path, default=None, help="Optional Django static destination for the generated PNG figures") args = parser.parse_args() df = build_dataset() models, train_idx, test_idx = evaluate(df) plot_figures(df, models, test_idx) lines = ['Synthetic classification of independent logistic-map sweeps', f'Independent sweep samples: {len(df)}', f'Train / test sweeps: {len(train_idx)} / {len(test_idx)}', 'Class 0 = hidden effective-r offset -0.045; class 1 = +0.045', 'Identical observed parameter grid for both classes', 'Only the training split is used to fit the classifiers and scalers.', 'pyGDIS is run independently on each complete sample sweep.', 'No claim of predictive benefit until independently replicated.',''] for name, m in models.items(): lines.append(f"{name}: balanced accuracy={m['balanced_accuracy']:.3f}, ROC AUC={m['auc']:.3f}") report = '\n'.join(lines) (OUT / 'metrics.txt').write_text(report + '\n', encoding='utf-8') print(report) print(f'Outputs: {OUT}') if args.publish_static_dir is not None: args.publish_static_dir.mkdir(parents=True, exist_ok=True) for png in OUT.glob('*.png'): shutil.copy2(png, args.publish_static_dir / png.name) print(f'Published PNG figures to: {args.publish_static_dir}') if __name__ == '__main__': main()