diff options
Diffstat (limited to 'examples/example_diffusion.py')
| -rw-r--r-- | examples/example_diffusion.py | 228 |
1 files changed, 228 insertions, 0 deletions
diff --git a/examples/example_diffusion.py b/examples/example_diffusion.py new file mode 100644 index 0000000..dc2d05e --- /dev/null +++ b/examples/example_diffusion.py @@ -0,0 +1,228 @@ +#!/usr/bin/env python3 +""" +Example demonstrating diffusion processing of transient events. + +This script shows how to use the new event_processor module to analyze +diffusion characteristics of detected events, similar to the approach +in heehun_diffusion_processing.py. +""" + +import numpy as np +import matplotlib.pyplot as plt + +# Import transivent functions +from transivent import ( + process_file, + get_final_events, + process_events_for_diffusion, + extract_event_waveforms, + calculate_msd_parallel, + calculate_acf, + fit_diffusion_linear, + plot_diffusion_comparison, + detect_events, + calculate_initial_background, + estimate_noise, + analyze_thresholds, +) + +# Configure logging +from transivent import configure_logging +configure_logging("INFO") + + +def create_synthetic_dataset( + n_events: int = 50, + event_duration: float = 100e-6, # 100 microseconds + sampling_interval: float = 1e-6, # 1 microsecond + total_duration: float = 0.1, # 100 milliseconds + diffusion_coeff: float = 1e-12, # m²/s + noise_level: float = 0.01, # Lower noise for better detection + signal_amplitude: float = 2.0, # Higher amplitude for better detection +): + """ + Create synthetic data with diffusion-like events. + + Parameters + ---------- + n_events : int + Number of events to generate + event_duration : float + Duration of each event in seconds + sampling_interval : float + Sampling interval in seconds + total_duration : float + Total duration of the signal in seconds + diffusion_coeff : float + Diffusion coefficient for event dynamics + noise_level : float + Noise level in the signal + signal_amplitude : float + Amplitude of the events + + Returns + ------- + t : np.ndarray + Time array + x : np.ndarray + Signal array + """ + n_points = int(total_duration / sampling_interval) + t = np.linspace(0, total_duration, n_points) + x = np.random.normal(0, noise_level, n_points) + + # Generate random event times + event_starts = np.random.uniform(0, total_duration - event_duration, n_events) + + for start_time in event_starts: + start_idx = int(start_time / sampling_interval) + end_idx = int((start_time + event_duration) / sampling_interval) + + if end_idx >= n_points: + end_idx = n_points - 1 + + # Create a diffusion-like signal for this event + event_length = end_idx - start_idx + event_t = np.linspace(0, event_duration, event_length) + + # Simulate Brownian motion + np.random.seed(int(start_time * 1e6)) # Seed for reproducibility + steps = np.random.normal(0, np.sqrt(2 * diffusion_coeff * sampling_interval), event_length) + event_signal = signal_amplitude * (1 + np.cumsum(steps)) + + # Add to main signal + x[start_idx:end_idx] += event_signal + + return t, x + + +def analyze_single_dataset(name: str, t: np.ndarray, x: np.ndarray, sampling_interval: float): + """ + Analyze a single dataset for diffusion characteristics. + + Parameters + ---------- + name : str + Name of the dataset + t : np.ndarray + Time array + x : np.ndarray + Signal array + sampling_interval : float + Sampling interval in seconds + + Returns + ------- + dict + Results containing events and diffusion analysis + """ + print(f"\n=== Analyzing {name} ===") + + # Step 1: Detect events using transivent + # First calculate background + smooth_n = int(10e-6 / sampling_interval) # 10 microsecond smoothing window + if smooth_n % 2 == 0: + smooth_n += 1 + + bg_initial = calculate_initial_background(t, x, smooth_n) + global_noise = estimate_noise(x, bg_initial) + + # Detect events + events, _ = detect_events( + t, x, bg_initial, + snr_threshold=3.0, + min_event_len=int(10e-6 / sampling_interval), + min_event_amp=5.0 * global_noise, + widen_frac=0.5, + global_noise=global_noise, + signal_polarity=1, # Positive events + ) + + print(f"Detected {len(events)} events") + + # Step 2: Process events for diffusion analysis + results = process_events_for_diffusion( + name=name, + sampling_interval=sampling_interval, + data_path="", + t=t, + x=x, + events=events, + max_lag=500, + n_jobs=-1, + ) + + print(f"Processed {results['event_count']} events for diffusion") + if results['event_count'] > 0: + print(f"Mean diffusion coefficient: {results['statistics']['mean_diffusion']:.3e} ± " + f"{results['statistics']['std_diffusion']:.3e} m²/s") + + return results + + +def main(): + """Main function demonstrating diffusion analysis.""" + print("Transivent Diffusion Processing Example") + print("=" * 40) + + # Create synthetic datasets with different diffusion characteristics + datasets = { + "Ferritin (slow diffusion)": { + "diffusion_coeff": 5e-13, # Slower diffusion + "signal_amplitude": 1.0, + "n_events": 40, + }, + "Catalase (medium diffusion)": { + "diffusion_coeff": 1e-12, # Medium diffusion + "signal_amplitude": 0.8, + "n_events": 35, + }, + "BSA (fast diffusion)": { + "diffusion_coeff": 2e-12, # Faster diffusion + "signal_amplitude": 0.6, + "n_events": 45, + }, + } + + sampling_interval = 1e-6 # 1 microsecond + + all_results = {} + + # Analyze each dataset + for name, params in datasets.items(): + # Generate synthetic data + t, x = create_synthetic_dataset( + n_events=params["n_events"], + diffusion_coeff=params["diffusion_coeff"], + signal_amplitude=params["signal_amplitude"], + sampling_interval=sampling_interval, + ) + + # Analyze for diffusion + results = analyze_single_dataset(name, t, x, sampling_interval) + all_results[name] = results + + # Create comparison plot + if any(len(results["diffusion_coeffs"]) > 0 for results in all_results.values()): + print("\n=== Creating comparison plot ===") + fig = plot_diffusion_comparison( + all_results, + show_ellipses=True, + colors=["#1f77b4", "gold", "green"], + markers=["o", "^", "s"], + ) + + # Save the plot + import os + os.makedirs("analysis_output", exist_ok=True) + fig.savefig("analysis_output/diffusion_comparison.png", dpi=150, bbox_inches="tight") + print("Saved plot to: analysis_output/diffusion_comparison.png") + + # Show plot + plt.show() + else: + print("\nNo events detected in any dataset. Check detection parameters.") + + +if __name__ == "__main__": + main()
\ No newline at end of file |