//! Shape optimization for aerodynamic design.

use aeroflow_shared::{
    AeroCoefficients, AirfoilGeometry, AnalysisType, FlowConditions, GeometryType, OperatorConfig,
    OptimizationConfig, OptimizationObjective, OptimizationResult, OptimizationStep,
    SimulationRequest,
};

use crate::geometry::CSTAirfoil;
use crate::AeroFlow;

/// Shape optimizer for airfoil design.
#[derive(Debug)]
pub struct ShapeOptimizer {
    /// Configuration.
    config: OptimizationConfig,
    /// RNG state.
    rng_state: u64,
}

impl ShapeOptimizer {
    /// Create a new shape optimizer.
    pub fn new(config: OptimizationConfig) -> Self {
        Self {
            config,
            rng_state: 42,
        }
    }

    /// Optimize airfoil shape.
    pub fn optimize(
        &mut self,
        aeroflow: &mut AeroFlow,
        initial_geometry: &AirfoilGeometry,
        conditions: &FlowConditions,
    ) -> OptimizationResult {
        // Convert to CST parameterization
        let mut cst = self.airfoil_to_cst(initial_geometry);
        let mut history = Vec::new();

        let mut best_objective = f32::MAX;
        let mut best_cst = cst.clone();
        let mut best_coeffs = AeroCoefficients::default();

        for iteration in 0..self.config.max_iterations {
            // Generate airfoil from CST
            let airfoil = cst.to_airfoil(50);

            // Evaluate aerodynamics
            let coeffs = self.evaluate_airfoil(aeroflow, &airfoil, conditions);

            // Compute objective
            let objective = self.compute_objective(&coeffs);
            let constraint_violation = self.check_constraints(&airfoil, &coeffs);

            // Record history
            history.push(OptimizationStep {
                iteration,
                objective,
                cl: coeffs.cl,
                cd: coeffs.cd,
                cm: coeffs.cm,
                constraint_violation,
            });

            // Update best if improved and feasible
            let penalized_objective = objective + 1000.0 * constraint_violation;
            if penalized_objective < best_objective {
                best_objective = penalized_objective;
                best_cst = cst.clone();
                best_coeffs = coeffs;
            }

            // Update design variables
            cst = self.gradient_step(&cst, aeroflow, conditions);
        }

        // Generate final geometry
        let optimized_geometry = best_cst.to_airfoil(100);

        OptimizationResult {
            optimized_geometry,
            objective_value: best_objective,
            coefficients: best_coeffs,
            history,
            num_evaluations: self.config.max_iterations,
            converged: best_objective < 1e10,
        }
    }

    /// Convert airfoil to CST parameterization.
    fn airfoil_to_cst(&self, airfoil: &AirfoilGeometry) -> CSTAirfoil {
        // Simplified conversion - in practice would fit CST coefficients
        let mut cst = CSTAirfoil::default();

        // Approximate based on thickness/camber
        let t = airfoil.max_thickness;
        let m = airfoil.max_camber;

        cst.upper_coeffs = vec![
            0.15 + t / 2.0 + m,
            0.25 + t / 3.0,
            0.20 + t / 4.0,
            0.10,
            0.05,
        ];
        cst.lower_coeffs = vec![-0.15 - t / 2.0 + m, -0.10 - t / 4.0, -0.08, -0.05, -0.02];

        cst
    }

    /// Evaluate aerodynamic coefficients for an airfoil.
    fn evaluate_airfoil(
        &self,
        aeroflow: &mut AeroFlow,
        airfoil: &AirfoilGeometry,
        conditions: &FlowConditions,
    ) -> AeroCoefficients {
        let request = SimulationRequest {
            geometry: GeometryType::Airfoil2D(airfoil.clone()),
            conditions: *conditions,
            operator: OperatorConfig::default(),
            analysis: AnalysisType::SinglePoint,
            compute_flow_field: false,
            export_results: false,
        };

        aeroflow.simulate(&request).coefficients
    }

    /// Compute objective function value.
    fn compute_objective(&self, coeffs: &AeroCoefficients) -> f32 {
        match &self.config.objective {
            OptimizationObjective::MinDragAtCl { target_cl } => {
                // Minimize drag with penalty for missing CL target
                let cl_penalty = 100.0 * (coeffs.cl - target_cl).powi(2);
                coeffs.cd + cl_penalty
            }
            OptimizationObjective::MaxLiftToDrag => {
                // Maximize L/D (minimize -L/D)
                -coeffs.cl / coeffs.cd.max(0.001)
            }
            OptimizationObjective::MaxLift => {
                // Maximize lift (minimize -CL)
                -coeffs.cl
            }
            OptimizationObjective::MinMomentVariation { cl_range } => {
                // Penalty for Cm variation
                let cl_in_range = coeffs.cl >= cl_range.0 && coeffs.cl <= cl_range.1;
                if cl_in_range {
                    coeffs.cm.abs()
                } else {
                    1000.0
                }
            }
            OptimizationObjective::Custom { weights } => {
                weights.drag_weight * coeffs.cd - weights.lift_weight * coeffs.cl
                    + weights.moment_weight * coeffs.cm.abs()
            }
        }
    }

    /// Check constraint violations.
    fn check_constraints(&self, airfoil: &AirfoilGeometry, coeffs: &AeroCoefficients) -> f32 {
        let mut violation = 0.0;

        // Thickness constraints
        if airfoil.max_thickness < self.config.constraints.min_thickness {
            violation += self.config.constraints.min_thickness - airfoil.max_thickness;
        }
        if airfoil.max_thickness > self.config.constraints.max_thickness {
            violation += airfoil.max_thickness - self.config.constraints.max_thickness;
        }

        // CL constraint
        if let Some(min_cl) = self.config.constraints.min_cl {
            if coeffs.cl < min_cl {
                violation += min_cl - coeffs.cl;
            }
        }

        // Cm constraint
        if let Some(max_cm) = self.config.constraints.max_cm_magnitude {
            if coeffs.cm.abs() > max_cm {
                violation += coeffs.cm.abs() - max_cm;
            }
        }

        violation
    }

    /// Perform gradient-based step.
    fn gradient_step(
        &mut self,
        cst: &CSTAirfoil,
        _aeroflow: &mut AeroFlow,
        _conditions: &FlowConditions,
    ) -> CSTAirfoil {
        // Simplified gradient descent with finite differences
        let mut new_cst = cst.clone();
        let lr = self.config.learning_rate;

        // Perturb upper surface coefficients
        for i in 0..new_cst.upper_coeffs.len() {
            new_cst.upper_coeffs[i] += lr * (self.random() as f32 - 0.5) * 0.02;
            // Clamp to reasonable range
            new_cst.upper_coeffs[i] = new_cst.upper_coeffs[i].clamp(-0.5, 0.5);
        }

        // Perturb lower surface coefficients
        for i in 0..new_cst.lower_coeffs.len() {
            new_cst.lower_coeffs[i] += lr * (self.random() as f32 - 0.5) * 0.02;
            new_cst.lower_coeffs[i] = new_cst.lower_coeffs[i].clamp(-0.5, 0.5);
        }

        new_cst
    }

    /// Random number generator.
    fn random(&mut self) -> f64 {
        self.rng_state = self
            .rng_state
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        (self.rng_state >> 11) as f64 / (1u64 << 53) as f64
    }
}

/// Genetic algorithm optimizer for global search.
#[derive(Debug)]
pub struct GeneticOptimizer {
    /// Population size.
    population_size: usize,
    /// Mutation rate.
    mutation_rate: f32,
    /// Crossover rate.
    crossover_rate: f32,
    /// RNG state.
    rng_state: u64,
}

impl GeneticOptimizer {
    /// Create a new genetic optimizer.
    pub fn new(population_size: usize) -> Self {
        Self {
            population_size,
            mutation_rate: 0.1,
            crossover_rate: 0.8,
            rng_state: 42,
        }
    }

    /// Initialize random population of CST airfoils.
    pub fn initialize_population(&mut self) -> Vec<CSTAirfoil> {
        (0..self.population_size)
            .map(|_| self.random_cst())
            .collect()
    }

    /// Generate random CST airfoil.
    fn random_cst(&mut self) -> CSTAirfoil {
        let mut cst = CSTAirfoil::default();

        for coeff in &mut cst.upper_coeffs {
            *coeff = (self.random() as f32 - 0.3) * 0.6;
        }
        for coeff in &mut cst.lower_coeffs {
            *coeff = (self.random() as f32 - 0.7) * 0.4;
        }

        cst
    }

    /// Tournament selection.
    pub fn tournament_select(&mut self, population: &[CSTAirfoil], fitness: &[f32]) -> CSTAirfoil {
        let tournament_size = 3;
        let mut best_idx = self.random_index(population.len());
        let mut best_fitness = fitness[best_idx];

        for _ in 1..tournament_size {
            let idx = self.random_index(population.len());
            if fitness[idx] < best_fitness {
                best_idx = idx;
                best_fitness = fitness[idx];
            }
        }

        population[best_idx].clone()
    }

    /// Crossover two parents.
    pub fn crossover(&mut self, parent1: &CSTAirfoil, parent2: &CSTAirfoil) -> CSTAirfoil {
        if self.random() as f32 > self.crossover_rate {
            return parent1.clone();
        }

        let mut child = parent1.clone();
        let crossover_point = self.random_index(child.upper_coeffs.len());

        for i in crossover_point..child.upper_coeffs.len() {
            child.upper_coeffs[i] = parent2.upper_coeffs[i];
            child.lower_coeffs[i] = parent2.lower_coeffs[i];
        }

        child
    }

    /// Mutate an individual.
    pub fn mutate(&mut self, individual: &mut CSTAirfoil) {
        for coeff in &mut individual.upper_coeffs {
            if (self.random() as f32) < self.mutation_rate {
                let delta = (self.random() as f32 - 0.5) * 0.1;
                *coeff += delta;
                *coeff = coeff.clamp(-0.5, 0.5);
            }
        }
        for coeff in &mut individual.lower_coeffs {
            if (self.random() as f32) < self.mutation_rate {
                let delta = (self.random() as f32 - 0.5) * 0.1;
                *coeff += delta;
                *coeff = coeff.clamp(-0.5, 0.5);
            }
        }
    }

    /// Random index.
    fn random_index(&mut self, max: usize) -> usize {
        (self.random() * max as f64) as usize % max
    }

    /// Random number generator.
    fn random(&mut self) -> f64 {
        self.rng_state = self
            .rng_state
            .wrapping_mul(6364136223846793005)
            .wrapping_add(1442695040888963407);
        (self.rng_state >> 11) as f64 / (1u64 << 53) as f64
    }
}

/// Adjoint-based gradient computation.
#[derive(Debug)]
pub struct AdjointGradient {
    /// Finite difference step size.
    epsilon: f32,
}

impl Default for AdjointGradient {
    fn default() -> Self {
        Self::new()
    }
}

impl AdjointGradient {
    /// Create a new adjoint gradient calculator.
    pub fn new() -> Self {
        Self { epsilon: 1e-4 }
    }

    /// Compute gradient via finite differences.
    pub fn compute_gradient(
        &self,
        aeroflow: &mut AeroFlow,
        cst: &CSTAirfoil,
        conditions: &FlowConditions,
        objective_fn: impl Fn(&AeroCoefficients) -> f32,
    ) -> CSTGradient {
        let mut upper_grad = vec![0.0; cst.upper_coeffs.len()];
        let mut lower_grad = vec![0.0; cst.lower_coeffs.len()];

        // Baseline evaluation
        let base_airfoil = cst.to_airfoil(50);
        let base_coeffs = self.evaluate(aeroflow, &base_airfoil, conditions);
        let base_obj = objective_fn(&base_coeffs);

        // Upper surface gradients
        for i in 0..cst.upper_coeffs.len() {
            let mut perturbed = cst.clone();
            perturbed.upper_coeffs[i] += self.epsilon;
            let perturbed_airfoil = perturbed.to_airfoil(50);
            let perturbed_coeffs = self.evaluate(aeroflow, &perturbed_airfoil, conditions);
            let perturbed_obj = objective_fn(&perturbed_coeffs);
            upper_grad[i] = (perturbed_obj - base_obj) / self.epsilon;
        }

        // Lower surface gradients
        for i in 0..cst.lower_coeffs.len() {
            let mut perturbed = cst.clone();
            perturbed.lower_coeffs[i] += self.epsilon;
            let perturbed_airfoil = perturbed.to_airfoil(50);
            let perturbed_coeffs = self.evaluate(aeroflow, &perturbed_airfoil, conditions);
            let perturbed_obj = objective_fn(&perturbed_coeffs);
            lower_grad[i] = (perturbed_obj - base_obj) / self.epsilon;
        }

        CSTGradient {
            upper: upper_grad,
            lower: lower_grad,
        }
    }

    /// Evaluate aerodynamic coefficients.
    fn evaluate(
        &self,
        aeroflow: &mut AeroFlow,
        airfoil: &AirfoilGeometry,
        conditions: &FlowConditions,
    ) -> AeroCoefficients {
        let request = SimulationRequest {
            geometry: GeometryType::Airfoil2D(airfoil.clone()),
            conditions: *conditions,
            operator: OperatorConfig::default(),
            analysis: AnalysisType::SinglePoint,
            compute_flow_field: false,
            export_results: false,
        };

        aeroflow.simulate(&request).coefficients
    }
}

/// Gradient of CST parameters.
#[derive(Debug, Clone)]
pub struct CSTGradient {
    /// Upper surface coefficient gradients.
    pub upper: Vec<f32>,
    /// Lower surface coefficient gradients.
    pub lower: Vec<f32>,
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_shape_optimizer_creation() {
        let config = OptimizationConfig::default();
        let optimizer = ShapeOptimizer::new(config);
        assert_eq!(optimizer.config.max_iterations, 100);
    }

    #[test]
    fn test_objective_max_ld() {
        let config = OptimizationConfig {
            objective: OptimizationObjective::MaxLiftToDrag,
            ..Default::default()
        };
        let optimizer = ShapeOptimizer::new(config);

        let coeffs = AeroCoefficients {
            cl: 0.5,
            cd: 0.01,
            ..Default::default()
        };

        let obj = optimizer.compute_objective(&coeffs);
        assert!((obj - (-50.0)).abs() < 0.01);
    }

    #[test]
    fn test_constraint_check() {
        let config = OptimizationConfig::default();
        let optimizer = ShapeOptimizer::new(config);

        let thin_airfoil = AirfoilGeometry {
            max_thickness: 0.05, // Below min_thickness of 0.06
            ..Default::default()
        };
        let coeffs = AeroCoefficients::default();

        let violation = optimizer.check_constraints(&thin_airfoil, &coeffs);
        assert!(violation > 0.0);
    }

    #[test]
    fn test_genetic_optimizer() {
        let mut ga = GeneticOptimizer::new(10);
        let population = ga.initialize_population();
        assert_eq!(population.len(), 10);
    }

    #[test]
    fn test_tournament_select() {
        let mut ga = GeneticOptimizer::new(10);
        let population = ga.initialize_population();
        let fitness: Vec<f32> = (0..10).map(|i| i as f32).collect();

        let selected = ga.tournament_select(&population, &fitness);
        // Should select from population
        assert_eq!(selected.upper_coeffs.len(), 5);
    }

    #[test]
    fn test_crossover() {
        let mut ga = GeneticOptimizer::new(10);
        ga.crossover_rate = 1.0; // Force crossover

        let parent1 = CSTAirfoil::default();
        let mut parent2 = CSTAirfoil::default();
        parent2.upper_coeffs = vec![0.5; 5];

        let child = ga.crossover(&parent1, &parent2);
        assert_eq!(child.upper_coeffs.len(), 5);
    }

    #[test]
    fn test_mutation() {
        let mut ga = GeneticOptimizer::new(10);
        ga.mutation_rate = 1.0; // Force mutation

        let mut individual = CSTAirfoil::default();
        let original = individual.upper_coeffs.clone();

        ga.mutate(&mut individual);

        // At least some coefficients should change
        let changed = individual
            .upper_coeffs
            .iter()
            .zip(original.iter())
            .any(|(a, b)| (a - b).abs() > 1e-6);
        assert!(changed);
    }

    #[test]
    fn test_adjoint_gradient() {
        let ag = AdjointGradient::new();
        assert!((ag.epsilon - 1e-4).abs() < 1e-6);
    }
}