//! # Liver Tumor Ablation Planning Example
//!
//! This example demonstrates using the Medical Digital Twin Platform to plan
//! radiofrequency ablation (RFA) treatment for a liver tumor.
//!
//! ## Workflow
//! 1. Create patient-specific liver geometry with tumor
//! 2. Create digital twin
//! 3. Run baseline simulation (no treatment)
//! 4. Plan RFA intervention
//! 5. Perform what-if analysis for different ablation parameters
//! 6. Select optimal treatment strategy
//!
//! ## Run
//! ```bash
//! cargo run --example liver_ablation_planning
//! ```

use rtx_digital_twin::{
    AblationProbe, DigitalTwin, InterventionType, OrganGeometry, TissueLabel, TissueType,
};

fn main() {
    println!("======================================");
    println!(" Liver Tumor Ablation Planning");
    println!("======================================\n");

    // Step 1: Create patient-specific geometry
    println!("[Step 1] Creating patient-specific liver geometry...");

    let shape = [60, 60, 60];
    let spacing = [1.0, 1.0, 1.0]; // 1mm voxel resolution
    let n = shape[0] * shape[1] * shape[2];
    let mut labels = vec![0u8; n]; // Start with air

    // Create liver parenchyma (realistic organ shape)
    for z in 10..50 {
        for y in 10..50 {
            for x in 10..50 {
                // Ellipsoidal liver shape
                let dx = (x as f32 - 30.0) / 20.0;
                let dy = (y as f32 - 30.0) / 20.0;
                let dz = (z as f32 - 30.0) / 20.0;
                let dist = dx * dx + dy * dy + dz * dz;

                if dist <= 1.0 {
                    let idx = z * shape[0] * shape[1] + y * shape[0] + x;
                    labels[idx] = TissueType::Liver.label();
                }
            }
        }
    }

    // Add tumor (3cm diameter at specific location)
    let tumor_center = [35, 30, 30];
    let tumor_radius = 8; // 8mm radius = 16mm diameter

    for z in 0..shape[2] {
        for y in 0..shape[1] {
            for x in 0..shape[0] {
                let dx = x as i32 - tumor_center[0] as i32;
                let dy = y as i32 - tumor_center[1] as i32;
                let dz = z as i32 - tumor_center[2] as i32;
                let dist_sq = dx * dx + dy * dy + dz * dz;

                if dist_sq <= tumor_radius * tumor_radius {
                    let idx = z * shape[0] * shape[1] + y * shape[0] + x;
                    if labels[idx] != 0 {
                        // Only replace liver, not air
                        labels[idx] = TissueType::Tumor.label();
                    }
                }
            }
        }
    }

    let geometry =
        OrganGeometry::from_labels(&labels, shape, spacing).expect("Failed to create geometry");

    // Print geometry statistics
    let histogram = geometry.tissue_histogram();
    let liver_voxels = histogram.get(&TissueType::Liver).unwrap_or(&0);
    let tumor_voxels = histogram.get(&TissueType::Tumor).unwrap_or(&0);

    println!("  - Grid size: {}x{}x{}", shape[0], shape[1], shape[2]);
    println!("  - Resolution: {:.1} mm", spacing[0]);
    println!(
        "  - Liver volume: ~{:.1} cm³",
        *liver_voxels as f32 / 1000.0
    );
    println!(
        "  - Tumor volume: ~{:.1} cm³",
        *tumor_voxels as f32 / 1000.0
    );
    println!();

    // Step 2: Create digital twin
    println!("[Step 2] Initializing digital twin...");
    let mut twin = DigitalTwin::new(geometry);

    let summary = twin.geometry_summary();
    println!("  - Total voxels: {}", summary.total_voxels);
    println!("  - Tissue voxels: {}", summary.tissue_voxels);
    println!(
        "  - Physical size: {:.1}x{:.1}x{:.1} mm",
        summary.dimensions[0], summary.dimensions[1], summary.dimensions[2]
    );
    println!();

    // Step 3: Baseline simulation
    println!("[Step 3] Running baseline simulation (no treatment)...");
    let baseline = twin
        .simulate_baseline()
        .expect("Baseline simulation failed");

    println!("  - Iterations: {}", baseline.iterations);
    println!("  - Max temperature: {:.1}°C", baseline.max_temperature);
    println!("  - Residual: {:.2e}", baseline.residual);
    println!();

    // Step 4: Plan RFA intervention
    println!("[Step 4] Planning RFA intervention...");
    println!(
        "  Target: Tumor center at ({}, {}, {})",
        tumor_center[0], tumor_center[1], tumor_center[2]
    );
    println!();

    // Step 5: What-if analysis for different power settings
    println!("[Step 5] What-if analysis: Comparing ablation strategies");
    println!("  Testing different power levels...\n");

    let power_levels = vec![30.0, 50.0, 70.0];
    let mut best_strategy = None;
    let mut best_score = f32::NEG_INFINITY;

    for (i, &power) in power_levels.iter().enumerate() {
        println!("  Scenario {}: {} W, 10 minutes", i + 1, power);

        // Reset between scenarios
        twin.reset_temperature();

        // Create probe
        let probe = AblationProbe::new(
            [
                tumor_center[0] as f32,
                tumor_center[1] as f32,
                tumor_center[2] as f32,
            ],
            power,
        )
        .with_type(InterventionType::RadiofrequencyAblation)
        .with_active_length(20.0) // 2cm active tip
        .with_diameter(3.0); // 3mm diameter

        // Run what-if analysis (10 minutes = 600 seconds)
        let result = twin
            .what_if(&probe, 600.0)
            .expect("What-if analysis failed");

        println!("    Max temp: {:.1}°C", result.max_temperature);
        println!(
            "    Damaged volume: {:.1} mm³ ({:.2} cm³)",
            result.total_damaged_volume,
            result.total_damaged_volume / 1000.0
        );
        println!("    Severe damage: {:.1} mm³", result.severe_damage_volume);
        println!(
            "    Safety OK: {}",
            if result.safety_margin_ok { "Yes" } else { "NO" }
        );
        println!("    Iterations: {}", result.iterations);

        // Simple scoring: want sufficient ablation but safe boundaries
        let tumor_volume_mm3 = *tumor_voxels as f32;
        let coverage = result.total_damaged_volume / tumor_volume_mm3;
        let safety_score = if result.safety_margin_ok { 1.0 } else { 0.0 };
        let score = coverage * 0.7 + safety_score * 0.3;

        println!("    Coverage: {:.1}%", coverage * 100.0);
        println!("    Score: {:.3}", score);
        println!();

        if score > best_score {
            best_score = score;
            best_strategy = Some((i + 1, power, result));
        }
    }

    // Step 6: Recommend optimal strategy
    println!("[Step 6] Treatment Recommendation");
    if let Some((scenario, power, result)) = best_strategy {
        println!("  Recommended: Scenario {} ({} W)", scenario, power);
        println!("  Expected outcomes:");
        println!("    - Max temperature: {:.1}°C", result.max_temperature);
        println!(
            "    - Ablation volume: {:.1} cm³",
            result.total_damaged_volume / 1000.0
        );
        println!(
            "    - Safety margin: {}",
            if result.safety_margin_ok {
                "Adequate"
            } else {
                "Review required"
            }
        );
        println!();

        // Generate clinical report
        println!("  Clinical Report:");
        println!("{}", result.report());
    }

    println!("======================================");
    println!(" Analysis Complete");
    println!("======================================");
}