📜 Curator: Visualizing Morphogenesis

.jules/curator.md

@@ -59,3 +59,11 @@
 **Gap:** The `freesurfer` module was a "Visual Void" with no explanation of the cortical reconstruction pipeline or usage.
 **Strategy:** Overhauled `freesurfer/mod.rs` with a Mermaid pipeline diagram and a runnable Quick Start example for cortical thickness calculation.
 **Outcome:** Users can now understand the MRI processing pipeline and use the tools for surface analysis.
+
+## 2026-02-02 - Visualizing Morphogenesis
+**Gap:** The `morphogenesis` module (Turing Patterns) is inherently visual but lacked diagrams explaining the numerical method (stencil) and had no runnable example, making it a "Visual Void".
+**Strategy:**
+- Added a Mermaid diagram to `morphogenesis.rs` visualizing the 1D Reaction-Diffusion stencil and sliding window optimization.
+- Added a "Quick Start" doc-test demonstrating system initialization and perturbation.
+- Created `examples/turing_patterns.rs` to generate ASCII art of the patterns.
+**Outcome:** Users can now understand the stencil logic and generate Turing patterns immediately.

math_explorer/examples/turing_patterns.rs

@@ -0,0 +1,104 @@
+//! # Turing Pattern Simulation Example
+//!
+//! This example simulates a 1D Reaction-Diffusion system using Schnakenberg kinetics.
+//! It demonstrates how to set up the system, perturb the initial state, and evolve it over time.
+//!
+//! Run with: `cargo run --example turing_patterns`
+
+use math_explorer::biology::morphogenesis::{SchnakenbergKinetics, TuringSystem};
+
+fn main() {
+    println!("🧪 Initializing Turing Pattern Simulation...");
+
+    // 1. Configuration
+    // Domain size: 60 points
+    // Diffusion coefficients: D_u = 1.0, D_v = 40.0 (Ratio 40)
+    let size = 60;
+    let d_u = 1.0;
+    let d_v = 40.0;
+    let dx = 1.0;
+
+    // Kinetics parameters
+    let a = 0.1;
+    let b = 0.9;
+    let kinetics = SchnakenbergKinetics::new(a, b);
+
+    // Calculate Homogeneous Steady State
+    // u* = a + b
+    // v* = b / (a + b)^2
+    let u_star = a + b;
+    let v_star = b / (u_star * u_star);
+
+    println!("   Steady State -> u*: {:.2}, v*: {:.2}", u_star, v_star);
+
+    let mut system =
+        TuringSystem::<SchnakenbergKinetics>::new_with_kinetics(size, d_u, d_v, dx, kinetics);
+
+    // 2. Initialize to Steady State + Perturbation
+    for i in 0..size {
+        system.u_mut()[i] = u_star;
+        system.v_mut()[i] = v_star;
+    }
+
+    // Perturb the center
+    let center = size / 2;
+    // Add random-ish noise (deterministic for this example)
+    for i in 0..size {
+        let noise = ((i as f64 * 0.1).sin()) * 0.01;
+        system.u_mut()[i] += noise;
+        system.v_mut()[i] += noise;
+    }
+    system.u_mut()[center] += 0.1; // Strong perturbation
+
+    println!("   Perturbation applied.");
+
+    // 3. Simulation Loop
+    let dt = 0.01;
+    let steps = 10000;
+
+    println!("   Simulating for {} steps (dt={})...", steps, dt);
+
+    for _ in 0..steps {
+        system.step(dt);
+    }
+
+    // 4. Output Analysis
+    println!("✅ Simulation Complete.");
+
+    // Simple ASCII visualization of the Activator (U)
+    println!("\n   Activator (U) Profile:");
+    print_ascii_profile(system.u());
+}
+
+fn print_ascii_profile(data: &[f64]) {
+    let max_val = data.iter().cloned().fold(f64::NEG_INFINITY, f64::max);
+    let min_val = data.iter().cloned().fold(f64::INFINITY, f64::min);
+
+    println!("   Min: {:.4}, Max: {:.4}", min_val, max_val);
+
+    if (max_val - min_val).abs() < 1e-6 {
+        println!("   [Flat Profile]");
+        return;
+    }
+
+    let range = max_val - min_val;
+    print!("   |");
+    for &val in data {
+        // Map value to 0..=9
+        let normalized = ((val - min_val) / range * 9.0).round() as usize;
+        let char = match normalized {
+            0 => ' ',
+            1 => '.',
+            2 => ':',
+            3 => '-',
+            4 => '=',
+            5 => '+',
+            6 => '*',
+            7 => '#',
+            8 => '%',
+            _ => '@',
+        };
+        print!("{}", char);
+    }
+    println!("|");
+}

math_explorer/src/biology/morphogenesis.rs

@@ -5,6 +5,61 @@
 //!
 //! The general equation is:
 //! $$ \frac{\partial \mathbf{u}}{\partial t} = D \nabla^2 \mathbf{u} + \mathbf{f}(\mathbf{u}) $$
+//!
+//! ## 🚀 Quick Start
+//!
+//! ```rust
+//! use math_explorer::biology::morphogenesis::{TuringSystem, SchnakenbergKinetics};
+//! use math_explorer::pure_math::analysis::ode::TimeStepper;
+//!
+//! // 1. Initialize System (Grid Size=100, Du=1.0, Dv=20.0, dx=1.0)
+//! // High ratio of diffusion coefficients (Dv >> Du) is required for patterns.
+//! let mut system = TuringSystem::<SchnakenbergKinetics>::new(100, 1.0, 20.0, 1.0);
+//!
+//! // 2. Perturb the center to break symmetry (crucial for pattern formation)
+//! system.u_mut()[50] += 1.0;
+//! system.v_mut()[50] += 1.0;
+//!
+//! // 3. Evolve for some time
+//! for _ in 0..100 {
+//!     system.step(0.01);
+//! }
+//!
+//! // 4. Inspect results
+//! let center_u = system.u()[50];
+//! println!("Center Activator Concentration: {:.4}", center_u);
+//! assert!(center_u != 0.0);
+//! ```
+//!
+//! ## ⚙️ Architecture: Stencil Optimization
+//!
+//! The solver uses a **Sliding Window** approach to compute the Laplacian $\nabla^2$ efficiently without
+//! unnecessary memory lookups. This allows the compiler to rotate registers for `u_prev`, `u_curr`, and `u_next`.
+//!
+//! ```mermaid
+//! graph LR
+//!     subgraph Memory[Grid Memory]
+//!     Cells[... u_{i-1}, u_{i}, u_{i+1} ...]
+//!     end
+//!
+//!     subgraph Registers[Sliding Window]
+//!     Prev[u_{i-1}]
+//!     Curr[u_{i}]
+//!     Next[u_{i+1}]
+//!     end
+//!
+//!     Cells -.->|Load| Next
+//!     Prev -->|Shift| Temp[Discard]
+//!     Curr -->|Shift| Prev
+//!     Next -->|Shift| Curr
+//!
+//!     Prev & Curr & Next -->|Finite Difference| Laplacian[∇²u]
+//!     Curr -->|Kinetics| Reaction[f(u,v)]
+//!
+//!     Laplacian & Reaction -->|Update| NewState[u_{i} + dt * du/dt]
+//!
+//!     style Curr fill:#f9f,stroke:#333,stroke-width:2px
+//! ```
 
 use crate::pure_math::analysis::ode::{OdeSystem, TimeStepper, VectorOperations};
 use std::ops::{Add, AddAssign, Mul, MulAssign};