KiThe 0.3.4

//! # Task Parser for Reactor BVP Module
//!
//! This module provides comprehensive parsing capabilities for chemical reactor boundary value problems (BVPs)
//! from structured text files. It bridges the gap between human-readable configuration files and the complex
//! data structures required by numerical BVP solvers.
//!
//! ## Purpose
//!
//! The module addresses the challenge of configuring complex reactor simulations through text files rather
//! than hardcoded parameters. It enables users to define complete reactor problems including:
//! - Chemical kinetics (reactions, Arrhenius parameters)
//! - Thermodynamic properties (heat capacity, thermal conductivity)
//! - Transport properties (diffusion coefficients)
//! - Boundary conditions and initial conditions
//! - Solver settings (tolerances, bounds, numerical methods)
//!
//! ## Main Methods
//!
//! - **`set_reactor_params_from_hashmap()`**: Core parsing method that extracts reactor parameters from DocumentMap
//! - **`parse_tolerance_and_bounds()`**: Specialized parser for solver tolerance and bounds configurations
//! - **`parse_file()`**: File I/O wrapper that creates DocumentParser from file path
//! - **`parse_parameters_with_exact_names()`**: High-level parsing orchestrator
//! - **`solve_from_map()`**: Complete workflow from parsed data to solved BVP
//! - **`solve_from_file()`**: One-shot method: file → parsing → solving
//!
//! ## Non-Obvious Code Features & Tips
//!
//! ### Variable Shadowing Prevention
//! The code carefully avoids variable name conflicts when parsing nested structures.
//! Example: `problem_name_val` instead of reusing `process_conditions` in if-let statements.
//!
//! ### DocumentMap Navigation Pattern
//! Uses consistent `.get().expect().clone().unwrap()[0].as_type()` pattern for extracting
//! values from the nested HashMap structure returned by the document parser.
//!
//! ### Groups Parsing Logic
//! The groups parsing implements a two-stage process: first check if groups are enabled,
//! then iterate through substances to build the nested HashMap structure.
//!
//! ### Vector Extraction Trick
//! For arrays like thermal_effects and reaction parameters, uses `.as_vector().unwrap().clone()`
//! to extract Vec<f64> directly from the parser's Value enum.
//!
//! ### Bounds Tuple Construction
//! Bounds parsing extracts 2-element vectors and converts them to tuples using array indexing:
//! `(bounds_vec[0], bounds_vec[1])` for efficient tuple creation.
//!
//! ### Error Handling Strategy
//! Uses `.expect()` for required parameters and `if let Some()` for optional ones,
//! providing clear error messages for missing required configuration.

use super::SimpleReactorBVP::{FastElemReact, SimpleReactorTask};
use crate::ReactorsBVP::reactor_BVP_utils::InitialConfig;
use crate::ReactorsBVP::reactor_BVP_utils::{BoundsConfig, ScalingConfig, ToleranceConfig};
use crate::ReactorsBVP::solver_backend::{
    ReactorBvpAotBuildPolicy, ReactorBvpAotBuildProfile, ReactorBvpAotCompilePreset,
    ReactorBvpAotCompiler, ReactorBvpAotExecutionPolicy, ReactorBvpMatrixBackend,
    ReactorBvpSolverConfig, ReactorBvpSymbolicBackend,
};
use crate::Utils::show_this_pic::show_image;
use RustedSciThe::command_interpreter::task_parser::{DocumentMap, DocumentParser, Value};
use RustedSciThe::numerical::BVP_Damp::NR_Damp_solver_damped::{AdaptiveGridConfig, SolverParams};
use RustedSciThe::numerical::BVP_Damp::grid_api::GridRefinementMethod;
use log::info;
use nalgebra::DMatrix;
use std::collections::HashMap;

fn missing_field(section: &str, field: &str) -> crate::ReactorsBVP::SimpleReactorBVP::ReactorError {
    crate::ReactorsBVP::SimpleReactorBVP::ReactorError::MissingData(format!(
        "Missing `{}` in `{}` section",
        field, section
    ))
}

fn invalid_numeric(
    section: &str,
    field: &str,
    value: f64,
) -> crate::ReactorsBVP::SimpleReactorBVP::ReactorError {
    crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidNumericValue(format!(
        "`{}.{} = {}` is not finite",
        section, field, value
    ))
}

fn required_section<'a>(
    task_hashmap: &'a DocumentMap,
    section: &str,
) -> Result<
    &'a HashMap<String, Option<Vec<Value>>>,
    crate::ReactorsBVP::SimpleReactorBVP::ReactorError,
> {
    task_hashmap
        .get(section)
        .ok_or_else(|| missing_field("document", section))
}

fn required_values<'a>(
    section: &'a HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<&'a [Value], crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    section
        .get(field)
        .ok_or_else(|| missing_field(section_name, field))?
        .as_ref()
        .map(|values| values.as_slice())
        .ok_or_else(|| missing_field(section_name, field))
}

fn required_f64(
    section: &HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<f64, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let value = required_values(section, section_name, field)?
        .first()
        .ok_or_else(|| missing_field(section_name, field))?
        .as_float()
        .ok_or_else(|| missing_field(section_name, field))?;
    if !value.is_finite() {
        return Err(invalid_numeric(section_name, field, value));
    }
    Ok(value)
}

fn required_usize(
    section: &HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<usize, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let value = required_values(section, section_name, field)?
        .first()
        .ok_or_else(|| missing_field(section_name, field))?;
    value_as_usize(value).ok_or_else(|| missing_field(section_name, field))
}

fn required_string(
    section: &HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<String, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    required_values(section, section_name, field)?
        .first()
        .ok_or_else(|| missing_field(section_name, field))?
        .as_string()
        .cloned()
        .ok_or_else(|| missing_field(section_name, field))
}

fn optional_string(
    section: &HashMap<String, Option<Vec<Value>>>,
    field: &str,
) -> Result<Option<String>, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    Ok(match section.get(field) {
        None => None,
        Some(values) => {
            let values = values
                .as_ref()
                .ok_or_else(|| missing_field("section", field))?;
            let value = values
                .first()
                .ok_or_else(|| missing_field("section", field))?;
            value
                .as_option_string()
                .cloned()
                .or_else(|| value.as_string().cloned())
        }
    })
}

// RustedSciThe's parser historically emits both Integer and Usize values.
// Keeping this conversion local avoids relying on one exact parser variant.
fn value_as_usize(value: &Value) -> Option<usize> {
    match value {
        Value::Usize(value) => Some(*value),
        Value::Integer(value) if *value >= 0 => Some(*value as usize),
        Value::Optional(Some(inner)) => value_as_usize(inner),
        _ => None,
    }
}

fn optional_usize(
    section: &HashMap<String, Option<Vec<Value>>>,
    field: &str,
) -> Result<Option<usize>, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    Ok(match section.get(field) {
        None => None,
        Some(values) => {
            let values = values
                .as_ref()
                .ok_or_else(|| missing_field("section", field))?;
            let value = values
                .first()
                .ok_or_else(|| missing_field("section", field))?;
            let parsed = value_as_usize(value)
                .filter(|value| *value > 0)
                .ok_or_else(|| {
                    crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration(
                        format!("`{}` must be a positive integer", field),
                    )
                })?;
            Some(parsed)
        }
    })
}

fn required_vector_f64(
    section: &HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<Vec<f64>, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let values = required_values(section, section_name, field)?;
    let vector = values
        .first()
        .ok_or_else(|| missing_field(section_name, field))?
        .as_vector()
        .ok_or_else(|| missing_field(section_name, field))?;
    if let Some(invalid) = vector.iter().copied().find(|value| !value.is_finite()) {
        return Err(invalid_numeric(section_name, field, invalid));
    }
    Ok(vector.clone())
}

fn required_bool(
    section: &HashMap<String, Option<Vec<Value>>>,
    section_name: &str,
    field: &str,
) -> Result<bool, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    required_values(section, section_name, field)?
        .first()
        .ok_or_else(|| missing_field(section_name, field))?
        .as_boolean()
        .ok_or_else(|| missing_field(section_name, field))
}

fn parse_tolerance_and_bounds_from_map(
    task_hashmap: &DocumentMap,
) -> Result<(ToleranceConfig, BoundsConfig), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let rel_tolerance = required_section(task_hashmap, "rel_tolerance")?;
    let tolerance_config = ToleranceConfig::new(
        required_f64(rel_tolerance, "rel_tolerance", "C")?,
        required_f64(rel_tolerance, "rel_tolerance", "J")?,
        required_f64(rel_tolerance, "rel_tolerance", "Teta")?,
        required_f64(rel_tolerance, "rel_tolerance", "q")?,
    );

    let bounds = required_section(task_hashmap, "bounds")?;
    let parse_pair =
        |field: &str| -> Result<(f64, f64), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
            let values = required_values(bounds, "bounds", field)?;
            if values.len() == 1 {
                let pair = values
                    .first()
                    .ok_or_else(|| missing_field("bounds", field))?
                    .as_vector()
                    .ok_or_else(|| missing_field("bounds", field))?;
                if pair.len() != 2 {
                    return Err(
                        crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration(
                            format!("`bounds.{}` must contain exactly two numbers", field),
                        ),
                    );
                }
                if let Some(invalid) = pair.iter().copied().find(|value| !value.is_finite()) {
                    return Err(invalid_numeric("bounds", field, invalid));
                }
                Ok((pair[0], pair[1]))
            } else if values.len() == 2 {
                let left = values[0]
                    .as_float()
                    .ok_or_else(|| missing_field("bounds", field))?;
                let right = values[1]
                    .as_float()
                    .ok_or_else(|| missing_field("bounds", field))?;
                if !left.is_finite() {
                    return Err(invalid_numeric("bounds", field, left));
                }
                if !right.is_finite() {
                    return Err(invalid_numeric("bounds", field, right));
                }
                Ok((left, right))
            } else {
                Err(
                    crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration(
                        format!(
                            "`bounds.{}` must contain one pair or two scalar values",
                            field
                        ),
                    ),
                )
            }
        };
    let bounds_config = BoundsConfig::new(
        parse_pair("C")?,
        parse_pair("J")?,
        parse_pair("Teta")?,
        parse_pair("q")?,
    );

    Ok((tolerance_config, bounds_config))
}

fn parse_basic_settings_from_map(
    task_hashmap: &DocumentMap,
) -> Result<(f64, f64, usize, String), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let process_conditions = required_section(task_hashmap, "process_conditions")?;
    let t0 = required_f64(process_conditions, "process_conditions", "t0")?;
    let t_end = required_f64(process_conditions, "process_conditions", "t_end")?;
    let n_steps = required_usize(process_conditions, "process_conditions", "n_steps")?;
    let arg = match process_conditions.get("arg") {
        Some(value) => {
            let values = value
                .as_ref()
                .ok_or_else(|| missing_field("process_conditions", "arg"))?;
            let arg = values
                .first()
                .ok_or_else(|| missing_field("process_conditions", "arg"))?
                .as_string()
                .cloned()
                .ok_or_else(|| missing_field("process_conditions", "arg"))?;
            Some(arg)
        }
        None => None,
    }
    .unwrap_or_else(|| "x".to_string());
    Ok((t0, t_end, n_steps, arg))
}

fn set_initial_guess_from_map_data(
    task_hashmap: &DocumentMap,
    n_steps: usize,
    unknown_count: usize,
) -> Result<DMatrix<f64>, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let inconfing = InitialConfig::new();
    let initial_guess: DMatrix<f64> = if let Some(guess_info) = task_hashmap.get("initial_guess") {
        let res = if let Some(universal) = guess_info.get("universal") {
            let values = universal
                .as_ref()
                .ok_or_else(|| missing_field("initial_guess", "universal"))?;
            let val = values
                .first()
                .ok_or_else(|| missing_field("initial_guess", "universal"))?
                .as_float()
                .ok_or_else(|| missing_field("initial_guess", "universal"))?;
            if !val.is_finite() {
                return Err(invalid_numeric("initial_guess", "universal", val));
            }
            inconfing.only_one_value_for_all_initial(val, n_steps, unknown_count)?
        } else {
            let C_val = required_f64(guess_info, "initial_guess", "C")?;
            let J_val = required_f64(guess_info, "initial_guess", "J")?;
            let Teta_val = required_f64(guess_info, "initial_guess", "Teta")?;
            let q_val = required_f64(guess_info, "initial_guess", "q")?;
            let map = HashMap::from([
                ("C".to_string(), C_val),
                ("J".to_string(), J_val),
                ("Teta".to_string(), Teta_val),
                ("q".to_string(), q_val),
            ]);
            inconfing.all_const_initial(map, n_steps)?
        };
        res
    } else {
        inconfing.only_one_value_for_all_initial(1e-2, n_steps, unknown_count)?
    };
    Ok(initial_guess)
}

#[derive(Debug)]
struct ParsedSolverSettings {
    solver_backend_config: ReactorBvpSolverConfig,
    scheme: String,
    strategy: String,
    strategy_params: Option<SolverParams>,
    linear_sys_method: Option<String>,
    method: String,
    abs_tolerance: f64,
    max_iterations: usize,
    loglevel: Option<String>,
    dont_save_log: bool,
}

fn parse_solver_settings_from_map(
    task_hashmap: &DocumentMap,
) -> Result<ParsedSolverSettings, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let solver_settings = required_section(task_hashmap, "solver_settings")?;
    let mut solver_backend_config =
        if let Some(generated_backend) = optional_string(solver_settings, "generated_backend")? {
            ReactorBvpSolverConfig::from_generated_backend_name(&generated_backend)
                .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?
        } else {
            ReactorBvpSolverConfig::default_lambdify()
        };
    if let Some(matrix_backend) = optional_string(solver_settings, "matrix_backend")? {
        solver_backend_config = solver_backend_config.with_matrix_backend(
            ReactorBvpMatrixBackend::parse(&matrix_backend).map_err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration,
            )?,
        );
    }
    if let Some(symbolic_backend) = optional_string(solver_settings, "symbolic_backend")? {
        solver_backend_config = solver_backend_config.with_symbolic_backend(
            ReactorBvpSymbolicBackend::parse(&symbolic_backend).map_err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration,
            )?,
        );
    }
    let aot_compiler = optional_string(solver_settings, "aot_compiler")?
        .or(optional_string(solver_settings, "aot_c_compiler")?);
    if let Some(compiler) = aot_compiler {
        let compiler = ReactorBvpAotCompiler::parse(&compiler)
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?;
        solver_backend_config =
            solver_backend_config.map_aot_config(|config| config.with_compiler(compiler));
    }
    if let Some(build_policy) = optional_string(solver_settings, "aot_build_policy")? {
        let build_policy = ReactorBvpAotBuildPolicy::parse(&build_policy)
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?;
        solver_backend_config =
            solver_backend_config.map_aot_config(|config| config.with_build_policy(build_policy));
    }
    if let Some(build_profile) = optional_string(solver_settings, "aot_build_profile")? {
        let build_profile = ReactorBvpAotBuildProfile::parse(&build_profile)
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?;
        solver_backend_config =
            solver_backend_config.map_aot_config(|config| config.with_build_profile(build_profile));
    }
    if let Some(compile_preset) = optional_string(solver_settings, "aot_compile_preset")? {
        let compile_preset = ReactorBvpAotCompilePreset::parse(&compile_preset)
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?;
        solver_backend_config = solver_backend_config
            .map_aot_config(|config| config.with_compile_preset(compile_preset));
    }
    if let Some(execution_policy) = optional_string(solver_settings, "aot_execution_policy")? {
        let execution_policy = ReactorBvpAotExecutionPolicy::parse(&execution_policy)
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration)?;
        solver_backend_config = solver_backend_config
            .map_aot_config(|config| config.with_execution_policy(execution_policy));
    }
    let target_chunks = optional_usize(solver_settings, "aot_target_chunks")?;
    let residual_chunks =
        optional_usize(solver_settings, "aot_residual_target_chunks")?.or(target_chunks);
    let jacobian_chunks =
        optional_usize(solver_settings, "aot_jacobian_target_chunks")?.or(target_chunks);
    if residual_chunks.is_some() || jacobian_chunks.is_some() {
        solver_backend_config = solver_backend_config
            .map_aot_config(|config| config.with_target_chunks(residual_chunks, jacobian_chunks));
    }

    let scheme = required_string(solver_settings, "solver_settings", "scheme")?;
    let strategy = required_string(solver_settings, "solver_settings", "strategy")?;
    let linear_sys_method = optional_string(solver_settings, "linear_sys_method")?;
    let method = required_string(solver_settings, "solver_settings", "method")?;
    let abs_tolerance = required_f64(solver_settings, "solver_settings", "abs_tolerance")?;
    let max_iterations = required_usize(solver_settings, "solver_settings", "max_iterations")?;
    let loglevel = optional_string(solver_settings, "loglevel")?;

    // Keep log suppression compatible with both historical spellings.
    let dont_save_log = if solver_settings.contains_key("dont_save_log") {
        required_bool(solver_settings, "solver_settings", "dont_save_log")?
    } else if solver_settings.contains_key("dont_save_logs") {
        required_bool(solver_settings, "solver_settings", "dont_save_logs")?
    } else {
        true
    };

    let strategy_params = if let Some(strategy_params) = task_hashmap.get("strategy_params") {
        let max_jac = strategy_params
            .get("max_jac")
            .and_then(|value| value.as_ref())
            .and_then(|values| values.first())
            .and_then(|value| value.as_option_usize());
        let max_damp_iter = strategy_params
            .get("max_damp_iter")
            .and_then(|value| value.as_ref())
            .and_then(|values| values.first())
            .and_then(|value| value.as_option_usize());
        let damp_factor = strategy_params
            .get("damp_factor")
            .and_then(|value| value.as_ref())
            .and_then(|values| values.first())
            .and_then(|value| value.as_option_float());

        let adaptive = if let Some(adaptive_strategy) = task_hashmap.get("adaptive_strategy") {
            let version = required_usize(adaptive_strategy, "adaptive_strategy", "version")?;
            let max_refinements =
                required_usize(adaptive_strategy, "adaptive_strategy", "max_refinements")?;
            let grid_refinement = required_section(task_hashmap, "grid_refinement")?;
            let (method_name, method_values) = grid_refinement
                .iter()
                .next()
                .ok_or_else(|| missing_field("grid_refinement", "method"))?;
            let method_params = method_values
                .as_ref()
                .ok_or_else(|| missing_field("grid_refinement", method_name))?
                .first()
                .ok_or_else(|| missing_field("grid_refinement", method_name))?
                .as_vector()
                .ok_or_else(|| missing_field("grid_refinement", method_name))?;
            if let Some(invalid) = method_params
                .iter()
                .copied()
                .find(|value| !value.is_finite())
            {
                return Err(invalid_numeric("grid_refinement", method_name, invalid));
            }
            let grid_method = match method_name.as_str() {
                "doubleoints" | "double_points" => GridRefinementMethod::DoublePoints,
                "easy" => GridRefinementMethod::Easy(method_params[0]),
                "grcarsmooke" | "grcar_smooke" => GridRefinementMethod::GrcarSmooke(
                    method_params[0],
                    method_params[1],
                    method_params[2],
                ),
                "pearson" => GridRefinementMethod::Pearson(method_params[0], method_params[1]),
                "twopnt" => GridRefinementMethod::TwoPoint(
                    method_params[0],
                    method_params[1],
                    method_params[2],
                ),
                _ => {
                    return Err(
                        crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration(
                            format!("Unknown grid refinement method `{}`", method_name),
                        ),
                    );
                }
            };
            Some(AdaptiveGridConfig {
                version,
                max_refinements,
                grid_method,
            })
        } else {
            None
        };

        Some(SolverParams {
            max_jac,
            max_damp_iter,
            damp_factor,
            adaptive,
        })
    } else {
        None
    };

    Ok(ParsedSolverSettings {
        solver_backend_config,
        scheme,
        strategy,
        strategy_params,
        linear_sys_method,
        method,
        abs_tolerance,
        max_iterations,
        loglevel,
        dont_save_log,
    })
}

/// Validate the dense solution matrix produced by the solver.
///
/// Solver backends may report a result without surfacing a structured failure.
/// This helper keeps the public contract strict by rejecting empty or non-finite
/// matrices before they are stored in reactor state.
fn validate_solver_output_matrix(
    solution: &DMatrix<f64>,
) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    if solution.nrows() == 0 || solution.ncols() == 0 {
        return Err(
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::CalculationError(
                "Solver returned an empty solution matrix".to_string(),
            ),
        );
    }
    if let Some((idx, value)) = solution
        .iter()
        .copied()
        .enumerate()
        .find(|(_, value)| !value.is_finite())
    {
        return Err(
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidNumericValue(format!(
                "Solver returned non-finite value {} at flat index {}",
                value, idx
            )),
        );
    }
    Ok(())
}

/// Build, configure, solve, and store a validated NRBVP snapshot.
///
/// The parser-facing solve methods share the same handoff logic; keeping it in
/// one place avoids drift between `solve_from_map` and `solve_from_parsed`.
fn solve_and_store_nrbvp(
    reactor: &mut SimpleReactorTask,
    task_map: &DocumentMap,
    tolerance_config: ToleranceConfig,
    bounds_config: BoundsConfig,
    t0: f64,
    t_end: f64,
    n_steps: usize,
    arg: String,
    initial_guess: DMatrix<f64>,
) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    let solver_settings = parse_solver_settings_from_map(task_map)?;

    let solver = &reactor.solver;
    let mut nr = solver.build_nrbvp_backend(
        crate::ReactorsBVP::SimpleReactorBVP::NrbvpHandoffConfig::new(
            initial_guess,
            t0,
            t_end,
            n_steps,
            solver_settings.scheme,
            solver_settings.strategy,
            solver_settings.strategy_params,
            solver_settings.linear_sys_method,
            solver_settings.method,
            solver_settings.abs_tolerance,
            Some(tolerance_config.to_full_tolerance_map(&reactor.kindata.substances)),
            solver_settings.max_iterations,
            Some(bounds_config.to_full_bounds_map(&reactor.kindata.substances)),
            solver_settings.loglevel,
            solver_settings.dont_save_log,
        )
        .with_solver_backend_config(solver_settings.solver_backend_config),
    )?;
    nr.before_solve_preprocessing();
    nr.solve();

    let solution = nr.get_result().ok_or_else(|| {
        crate::ReactorsBVP::SimpleReactorBVP::ReactorError::CalculationError(
            "Solver finished without producing a solution matrix".to_string(),
        )
    })?;
    validate_solver_output_matrix(&solution)?;
    if nr.x_mesh.is_empty() {
        return Err(
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::CalculationError(
                "Solver finished without producing an x mesh".to_string(),
            ),
        );
    }
    if nr.x_mesh.iter().any(|value| !value.is_finite()) {
        return Err(
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidNumericValue(
                "Solver x mesh contains non-finite values".to_string(),
            ),
        );
    }

    reactor.solver.x_mesh = Some(nr.x_mesh);
    reactor.solver.solution = Some(solution);
    reactor.solver.arg_name = arg;
    set_postprocessing_from_map(reactor, task_map)?;
    Ok(())
}

/// Apply post-processing options from a parsed document map.
///
/// This keeps the post-processing stage independent from `DocumentParser`
/// so the solver path can work from either parsed or already-materialized
/// configuration data without reparsing or mutating parser state.
fn set_postprocessing_from_map(
    reactor: &mut SimpleReactorTask,
    task_hashmap: &DocumentMap,
) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    reactor.check_balances()?;
    let solver_settings = required_section(task_hashmap, "postprocessing")?;

    // Flag to create plot via Rust native crate.
    let plot_flag = match solver_settings.get("plot") {
        Some(_) => required_bool(solver_settings, "postprocessing", "plot")?,
        None => false,
    };
    // Flag to create plot via GNU plot library.
    let gnuplot_flag = match solver_settings.get("gnuplot") {
        Some(_) => required_bool(solver_settings, "postprocessing", "gnuplot")?,
        None => false,
    };
    // Flag to save solution to txt.
    let save_flag = match solver_settings.get("save") {
        Some(_) => required_bool(solver_settings, "postprocessing", "save")?,
        None => false,
    };
    let save_to_csv = match solver_settings.get("save_to_csv") {
        Some(_) => required_bool(solver_settings, "postprocessing", "save_to_csv")?,
        None => false,
    };

    let name = optional_string(solver_settings, "filename")?;
    // Return from dimensionless to dimensioned unknowns.
    let return_to_dimension = match solver_settings.get("return_to_dimension") {
        Some(_) => required_bool(solver_settings, "postprocessing", "return_to_dimension")?,
        None => true,
    };
    let no_plots_in_terminal = match solver_settings.get("no_plots_in_terminal") {
        Some(_) => required_bool(solver_settings, "postprocessing", "no_plots_in_terminal")?,
        None => false,
    };

    if return_to_dimension {
        reactor.postprocessing()?;
    }
    if plot_flag {
        reactor.plot()?;
        let _ = show_image("Teta");
    }
    if gnuplot_flag {
        reactor.gnuplot()?;
        let _ = show_image("Teta");
    }
    if !no_plots_in_terminal {
        reactor.plot_in_terminal()?;
    }
    if save_flag {
        reactor.save_to_file(name.clone())?
    }
    if save_to_csv {
        reactor.save_to_csv(name)?;
    }
    reactor.estimate_values()?;
    Ok(())
}
impl SimpleReactorTask {
    /// Parses reactor parameters from a DocumentMap and populates the SimpleReactorTask structure.
    ///
    /// This is the core parsing method that extracts all reactor configuration from the parsed
    /// document structure. It handles reactions, process conditions, boundary conditions,
    /// diffusion coefficients, and chemical groups.
    ///
    /// # Arguments
    /// * `task_hashmap` - The parsed document structure containing all configuration sections
    ///
    /// # Sections Parsed
    /// - `reactions`: Arrhenius parameters (A, n, E, Q) for each reaction
    /// - `process_conditions`: Physical properties (Tm, P, Cp, Lambda, m, etc.)
    /// - `boundary_condition`: Initial/boundary values for all variables
    /// - `diffusion_coefficients`: Transport properties for each substance
    /// - Chemical groups: Atomic composition for custom substance definitions
    pub fn set_reactor_params_from_hashmap(
        &mut self,
        task_hashmap: &DocumentMap,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let reactions = required_section(task_hashmap, "reactions")?;
        let mut vec_of_struct: Vec<FastElemReact> = Vec::with_capacity(reactions.len());
        for (key, value) in reactions.iter() {
            let values = value
                .as_ref()
                .ok_or_else(|| missing_field("reactions", key))?;
            let rates = values
                .first()
                .ok_or_else(|| missing_field("reactions", key))?
                .as_vector()
                .ok_or_else(|| missing_field("reactions", key))?;
            if rates.len() != 4 {
                return Err(
                    crate::ReactorsBVP::SimpleReactorBVP::ReactorError::InvalidConfiguration(
                        format!(
                            "Reaction `{}` must store exactly 4 Arrhenius parameters",
                            key
                        ),
                    ),
                );
            }
            if let Some(invalid) = rates.iter().copied().find(|value| !value.is_finite()) {
                return Err(invalid_numeric("reactions", key, invalid));
            }
            vec_of_struct.push(FastElemReact {
                eq: key.to_owned(),
                A: rates[0],
                n: rates[1],
                E: rates[2],
                Q: rates[3],
            });
        }
        self.fast_react_set(vec_of_struct)?;

        let process_conditions = required_section(task_hashmap, "process_conditions")?;
        self.problem_name = optional_string(process_conditions, "problem_name")?;
        self.problem_description = optional_string(process_conditions, "problem_description")?;

        let substances_values =
            required_values(process_conditions, "process_conditions", "substances")?;
        let mut substances: Vec<String> = Vec::with_capacity(substances_values.len());
        for val in substances_values.iter() {
            let substance = val
                .as_string()
                .cloned()
                .ok_or_else(|| missing_field("process_conditions", "substances"))?;
            substances.push(substance);
        }
        self.kindata.substances = substances;

        self.Tm = required_f64(process_conditions, "process_conditions", "Tm")?;
        let L = required_f64(process_conditions, "process_conditions", "L")?;
        let dT = required_f64(process_conditions, "process_conditions", "dT")?;
        let T_scale = required_f64(process_conditions, "process_conditions", "T_scale")?;
        self.scaling = ScalingConfig::new(dT, L, T_scale);
        self.P = required_f64(process_conditions, "process_conditions", "P")?;
        self.Cp = required_f64(process_conditions, "process_conditions", "Cp")?;
        self.Lambda = required_f64(process_conditions, "process_conditions", "Lambda")?;
        self.m = required_f64(process_conditions, "process_conditions", "m")?;
        if let Some(mass) = process_conditions.get("M") {
            let values = mass
                .as_ref()
                .ok_or_else(|| missing_field("process_conditions", "M"))?;
            let value = values
                .first()
                .ok_or_else(|| missing_field("process_conditions", "M"))?
                .as_float()
                .ok_or_else(|| missing_field("process_conditions", "M"))?;
            if !value.is_finite() {
                return Err(invalid_numeric("process_conditions", "M", value));
            }
            self.M = value;
        }

        self.thermal_effects =
            required_vector_f64(process_conditions, "process_conditions", "thermal_effects")?;

        let boundary_condition_section = required_section(task_hashmap, "boundary_condition")?;
        let mut boundary_condition: HashMap<String, f64> =
            HashMap::with_capacity(boundary_condition_section.len());
        for (key, value) in boundary_condition_section.iter() {
            let values = value
                .as_ref()
                .ok_or_else(|| missing_field("boundary_condition", key))?;
            let value = values
                .first()
                .ok_or_else(|| missing_field("boundary_condition", key))?
                .as_float()
                .ok_or_else(|| missing_field("boundary_condition", key))?;
            if !value.is_finite() {
                return Err(invalid_numeric("boundary_condition", key, value));
            }
            boundary_condition.insert(key.to_owned(), value);
        }
        self.boundary_condition = boundary_condition;

        let diffusion_coefficients_section =
            required_section(task_hashmap, "diffusion_coefficients")?;
        let mut diffusion_coefficients: HashMap<String, f64> =
            HashMap::with_capacity(diffusion_coefficients_section.len());
        for (key, value) in diffusion_coefficients_section.iter() {
            let values = value
                .as_ref()
                .ok_or_else(|| missing_field("diffusion_coefficients", key))?;
            let value = values
                .first()
                .ok_or_else(|| missing_field("diffusion_coefficients", key))?
                .as_float()
                .ok_or_else(|| missing_field("diffusion_coefficients", key))?;
            if !value.is_finite() {
                return Err(invalid_numeric("diffusion_coefficients", key, value));
            }
            diffusion_coefficients.insert(key.to_owned(), value);
        }
        self.Diffusion = diffusion_coefficients;

        if let Some(groups_val) = process_conditions.get("groups") {
            let values = groups_val
                .as_ref()
                .ok_or_else(|| missing_field("process_conditions", "groups"))?;
            let groups_enabled = values
                .first()
                .ok_or_else(|| missing_field("process_conditions", "groups"))?
                .as_boolean()
                .ok_or_else(|| missing_field("process_conditions", "groups"))?;
            if groups_enabled {
                let mut groups: HashMap<String, HashMap<String, usize>> =
                    HashMap::with_capacity(self.kindata.substances.len());
                for sub_i in self.kindata.substances.iter() {
                    if let Some(sub_i_groups) = task_hashmap.get(sub_i.as_str()) {
                        let mut sub_i_groups_hashmap: HashMap<String, usize> =
                            HashMap::with_capacity(sub_i_groups.len());
                        for (key, value) in sub_i_groups.iter() {
                            let values = value.as_ref().ok_or_else(|| missing_field(sub_i, key))?;
                            let value = values
                                .first()
                                .ok_or_else(|| missing_field(sub_i, key))?
                                .as_usize()
                                .ok_or_else(|| missing_field(sub_i, key))?;
                            sub_i_groups_hashmap.insert(key.to_owned(), value);
                        }
                        groups.insert(sub_i.clone(), sub_i_groups_hashmap);
                    }
                }
                self.kindata.groups = Some(groups);
            }
        }

        Ok(())
    }
    /// Parses solver tolerance and bounds configurations from the document.
    ///
    /// Extracts numerical solver settings including relative tolerances and variable bounds
    /// for the four main variable types: C (concentrations), J (fluxes), Teta (temperature), q (heat flux).
    ///
    /// # Arguments
    /// * `parser` - DocumentParser containing the parsed configuration
    ///
    /// # Returns
    /// * `Ok((ToleranceConfig, BoundsConfig))` - Parsed tolerance and bounds configurations
    /// * `Err(String)` - Error message if required sections are missing
    ///
    /// # Expected Document Structure
    /// ```text
    /// rel_tolerance
    /// C: 1e-5
    /// J: 1e-5
    /// Teta: 1e-5
    /// q: 1e-5
    /// bounds
    /// C: -10.0, 10.0
    /// J: -1e20, 1e20
    /// Teta: -100.0, 100.0
    /// q: -1e20, 1e20
    /// ```
    fn parse_tolerance_and_bounds(
        &mut self,
        parser: &mut DocumentParser,
    ) -> Result<(ToleranceConfig, BoundsConfig), crate::ReactorsBVP::SimpleReactorBVP::ReactorError>
    {
        let task_hashmap = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "document must be parsed before tolerance extraction".to_string(),
            )
        })?;
        parse_tolerance_and_bounds_from_map(task_hashmap)
    }

    /// Legacy spelling kept for compatibility with older call sites.
    fn parse_toleranse_and_bounds(
        &mut self,
        parser: &mut DocumentParser,
    ) -> Result<(ToleranceConfig, BoundsConfig), crate::ReactorsBVP::SimpleReactorBVP::ReactorError>
    {
        self.parse_tolerance_and_bounds(parser)
    }

    /// Extracts basic solver settings (time domain and grid size) from the parsed document.
    ///
    /// # Arguments
    /// * `parser` - DocumentParser containing the configuration
    ///
    /// # Returns
    /// * `(t0, t_end, n_steps)` - Start time, end time, and number of grid steps
    pub fn parse_basic_settings(
        &mut self,
        parser: DocumentParser,
    ) -> Result<(f64, f64, usize, String), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let task_hashmap = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "document must be parsed before settings extraction".to_string(),
            )
        })?;
        parse_basic_settings_from_map(task_hashmap)
    }

    pub fn set_postprocessing_from_hashmap(
        &mut self,
        parser: &mut DocumentParser,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let result = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "document must be parsed before postprocessing".to_string(),
            )
        })?;
        set_postprocessing_from_map(self, result)
    }

    /// Legacy spelling kept for compatibility with older call sites.
    pub fn set_postpocessing_from_hashmap(
        &mut self,
        parser: &mut DocumentParser,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        self.set_postprocessing_from_hashmap(parser)
    }

    pub fn set_initial_guess_from_map(
        &mut self,
        parser: DocumentParser,
        n_steps: usize,
    ) -> Result<DMatrix<f64>, crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let result = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "document must be parsed before initial guess extraction".to_string(),
            )
        })?;
        set_initial_guess_from_map_data(result, n_steps, self.solver.unknowns.len())
    }
    /// Complete workflow: parses configuration, sets up BVP, configures solver, and solves the problem.
    ///
    /// This is a high-level method that orchestrates the entire process from parsed document
    /// to solved BVP. It handles all intermediate steps including BVP setup, solver configuration,
    /// tolerance/bounds parsing, and postprocessing.
    ///
    /// # Arguments
    /// * `parser` - DocumentParser containing the complete problem configuration
    ///
    /// # Process Flow
    /// 1. Parse reactor parameters
    /// 2. Setup BVP equations and boundary conditions
    /// 3. Configure NRBVP solver with parsed settings
    /// 4. Parse tolerances and bounds
    /// 5. Set initial guess and solve
    /// 6. Apply postprocessing (plotting, data export)
    pub fn solve_from_map(
        &mut self,
        mut parser: DocumentParser,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        self.parse_parameters_with_exact_names(&mut parser)?;
        let task_map = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "document must be parsed before solver setup".to_string(),
            )
        })?;
        let (tolerance_config, bounds_config) = parse_tolerance_and_bounds_from_map(task_map)?;
        let (t0, t_end, n_steps, arg) = parse_basic_settings_from_map(task_map)?;
        self.setup_bvp()?;
        if self.solver.unknowns.is_empty() {
            return Err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::MissingData(
                    "Solver unknowns were not initialized".to_string(),
                ),
            );
        }
        if self.solver.eq_system.is_empty() {
            return Err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::MissingData(
                    "Solver equation system was not initialized".to_string(),
                ),
            );
        }
        let initial_guess =
            set_initial_guess_from_map_data(task_map, n_steps, self.solver.unknowns.len())?;
        solve_and_store_nrbvp(
            self,
            task_map,
            tolerance_config,
            bounds_config,
            t0,
            t_end,
            n_steps,
            arg,
            initial_guess,
        )
    }

    /// One-shot method: reads configuration file, parses it, and solves the BVP problem.
    ///
    /// Convenience method that combines file reading, parsing, and solving in a single call.
    /// Ideal for batch processing or when the entire workflow can be automated.
    ///
    /// # Arguments
    /// * `path` - Path to the configuration file
    pub fn solve_from_file(
        &mut self,
        path: std::path::PathBuf,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let mut parser = self
            .parse_file(Some(path))
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError)?;
        parser
            .parse_document()
            .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError)?;
        self.solve_from_map(parser)
    }

    /// Creates a DocumentParser from a file path.
    ///
    /// Simple wrapper around DocumentParser creation that handles file I/O.
    ///
    /// # Arguments
    /// * `path` - Optional path to configuration file
    ///
    /// # Returns
    /// * `Ok(DocumentParser)` - Configured parser ready for document parsing
    /// * `Err(String)` - File I/O error message
    pub fn parse_file(
        &mut self,
        path: Option<std::path::PathBuf>,
    ) -> Result<DocumentParser, String> {
        let mut parser = DocumentParser::new(String::new());
        parser.setting_from_file(path)?;
        Ok(parser)
    }

    /// High-level parsing orchestrator that processes the document and populates reactor parameters.
    ///
    /// Coordinates the document parsing process and delegates to `set_reactor_params_from_hashmap`
    /// for actual parameter extraction.
    ///
    /// # Arguments
    /// * `parser` - DocumentParser with loaded configuration
    ///
    /// # Returns
    /// * `Ok(())` - Successful parsing
    /// * `Err(String)` - Parsing error message
    pub fn parse_parameters_with_exact_names(
        &mut self,
        parser: &mut DocumentParser,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        // Keep the helper idempotent: callers may hand us a parser that is already parsed.
        // Re-parsing can mutate the internal document shape in some backends, so we only
        // parse when the result cache is still empty.
        if parser.get_result().is_none() {
            parser
                .parse_document()
                .map_err(crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError)?;
        }
        let result = parser.get_result().ok_or_else(|| {
            crate::ReactorsBVP::SimpleReactorBVP::ReactorError::ParseError(
                "No result after parsing".to_string(),
            )
        })?;
        self.set_reactor_params_from_hashmap(result)?;
        Ok(())
    }

    pub fn solve_from_parsed(
        &mut self,
        map: DocumentMap,
    ) -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
        let task_map = map;
        self.set_reactor_params_from_hashmap(&task_map)?;
        let (tolerance_config, bounds_config) = parse_tolerance_and_bounds_from_map(&task_map)?;
        let (t0, t_end, n_steps, arg) = parse_basic_settings_from_map(&task_map)?;
        self.setup_bvp()?;
        if self.solver.unknowns.is_empty() {
            return Err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::MissingData(
                    "Solver unknowns were not initialized".to_string(),
                ),
            );
        }
        if self.solver.eq_system.is_empty() {
            return Err(
                crate::ReactorsBVP::SimpleReactorBVP::ReactorError::MissingData(
                    "Solver equation system was not initialized".to_string(),
                ),
            );
        }
        let initial_guess =
            set_initial_guess_from_map_data(&task_map, n_steps, self.solver.unknowns.len())?;
        solve_and_store_nrbvp(
            self,
            &task_map,
            tolerance_config,
            bounds_config,
            t0,
            t_end,
            n_steps,
            arg,
            initial_guess,
        )
    }
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
/// Creates a template configuration file for reactor BVP tasks.
///
/// Generates a commented template file that users can fill with their specific values.
/// The template includes all required sections with example values and explanatory comments.
///
pub const SIMPLE_BVP_TEMPLATE: &'static str = r#"
        initial_guess
        universal:1e-2
        process_conditions
        problem_name: Some(some_name)
        problem_description: Some("some decription")
        substances: A, B
        t0: 0.0
        t_end: 1.0
        n_steps: 25
        arg:x
        Tm: 1500.0
        L: 5e-4
        dT: 600.0
        T_scale: 600.0
        P: 1e6
        Cp: 1464.4
        Lambda: 0.07
        m: 0.0043
        M: 0.0342
        thermal_effects: [102000.0]
        groups:true
        boundary_condition
        A: 0.999
        B: 0.001
        T: 800.0
        diffusion_coefficients
        A: 0.000009296
        B: 0.000009296
        group1
        H: 4
        N: 8
        C: 8
        O: 8
        group2
        H: 6
        C: 1
        O: 1
        reactions
        A=>10B: [130000.0, 0.0, 20920.0, 102000.0]
        solver_settings
        scheme: forward
        method: Sparse
        strategy: Damped
        linear_sys_method: None
        abs_tolerance: 1e-7
        max_iterations: 100
        loglevel: Some(info)
        dont_save_logs: true
        bounds
        C: -10.0, 10.0
        J:  -1e20, 1e20
        Teta:-100.0, 100.0
        q: -1e20, 1e20
        rel_tolerance
        C: 1e-7
        J: 1e-7
        Teta: 1e-7
        q:  1e-7
        strategy_params
        max_jac: Some(3)
        max_damp_iter: Some(10)
        damp_factor: Some(0.5)
        # Adaptive grid refinement settings (optional)
        adaptive_strategy
        # Refinement version
        version: 1
        # Maximum refinement iterations
        max_refinements: 3
                
        #Grid refinement method and parameters
        grid_refinement
        // Available methods:
        // double_points: []
        // easy: [parameter]
        // grcar_smooke: [param1, param2, param3] (legacy token: grcarsmooke)
        // pearson: [param1, param2]
        // twopnt: [param1, param2, param3]
        grcar_smooke: [0.05, 0.05, 1.25]
        postprocessing
        gnuplot:true
        save_to_csv:false
        filename: meow
        gui_plot: true
        "#;
pub fn create_template() -> Result<(), crate::ReactorsBVP::SimpleReactorBVP::ReactorError> {
    use std::fs::File;
    use std::io::Write;

    let template_content = r#"
# Reactor BVP Configuration Template
// This template provides a starting point for configuring reactor BVP tasks.
// Users should fill in the values below for their specific problem.
// The template includes all required sections with example values and explanatory comments.
//
// some fields are filled with Some(value) it means the value in this field
// will be used in the task, otherwise the default value will be used
// if the field is None, the default value will be used
//
# Fill in the values below for your specific problem
# initial guess - if not set, default is 1e-2 for all variables and mesh
# if set universal:some_value - for all variables and mesh will be set this value
# if set C: some_value, J: some_value, Teta: some_value, q: some_value - for all 
# corresponing variables will be set this values
initial_guess:
universal: 1e-2
# Process conditions - main problem parameters
process_conditions
# Optional problem identification
problem_name: Some(YourProblemName)
problem_description: Some(YourProblemDescription)
# List of chemical substances (comma-separated)
substances: Substance1, Substance2
# Start position (real length is set in L field, so 
# t0 and t_end are dimensionless and reasonable to leave 
# them 0.0 and 1.0 )
t0: 0.0
# End position 
t_end: 1.0
# Number of grid points - the more you set 
# the more accurate the solution will be, but 
# the more time it will take to solve the problem
n_steps: 30
# Independent variable name (x, t, etc.)
arg: x
# Maximum temperature in the front of reaction [K]
Tm: 1500.0
# Characteristic length [m] - should be chosen in such 
# way that all reactions are ended at the end of the distance
L: 9e-4
# Temperature difference [K]
dT: 600.0
# Temperature scaling [K]
T_scale: 600.0
# Pressure [Pa]
P: 1e6
# Heat capacity [J/kg/K]
Cp: 1464.4
# Thermal conductivity [W/m/K]
Lambda: 0.07
# Mass velocity [kg/m2*s]
m: 0.0043
# Molar mass [kg/mol]
M: 0.0342
# Thermal effects for each reaction [J/mol]
thermal_effects: [102000.0]
# Enable custom atomic groups (true/false)
groups: true

# Boundary conditions - initial values for all variables
boundary_condition
# Initial concentration/mole fraction
Substance1: 0.999
# Initial concentration/mole fraction
Substance2: 0.001
# Initial temperature [K]
T: 800.0

# Diffusion coefficients for each substance [m²/s]
diffusion_coefficients
Substance1: 0.000009296
Substance2: 0.000009296

# Atomic composition (only if groups: true)
# Define atomic composition for each substance
Substance1
# Number of hydrogen atoms
H: 4
# Number of nitrogen atoms
N: 8
# Number of carbon atoms
C: 8
# Number of oxygen atoms
O: 8

Substance2
H: 6
C: 1
O: 1

# Chemical reactions with Arrhenius parameters
# Format: Reactant=>Product: [A, n, E, Q]
# A: pre-exponential factor, n: temperature exponent
# E: activation energy [J/mol/K],
# Q: heat of reaction [J/mol]
reactions
Substance1=>Substance2: [130000.0, 0.0, 20920.0, 102000.0]

# Numerical solver settings
solver_settings
# Discretization scheme
scheme: forward
# Solution method
method: Sparse
# Convergence strategy
strategy: Damped
# Linear system solver
linear_sys_method: None
# Absolute tolerance
abs_tolerance: 1e-5
# Maximum iterations
max_iterations: 100
# Logging level
loglevel: Some(info)
# Disable log saving
dont_save_logs: true

# Variable bounds for solver
bounds
# Concentration bounds
C: -10.0, 10.0
# Flux bounds
J: -1e20, 1e20
# Temperature bounds
Teta: -100.0, 100.0
# Heat flux bounds
q: -1e20, 1e20

# Relative tolerances for each variable type
rel_tolerance
# Concentration tolerance
C: 1e-5
# Flux tolerance
J: 1e-5
# Temperature tolerance
Teta: 1e-5
# Heat flux tolerance
q: 1e-5

# Advanced solver strategy parameters
strategy_params
# Maximum Jacobian updates
max_jac: Some(3)
# Maximum damping iterations
max_damp_iter: Some(10)
# Damping factor
damp_factor: Some(0.5)

# Adaptive grid refinement settings (optional)
adaptive_strategy
# Refinement version
version: 1
# Maximum refinement iterations
max_refinements: 3
        
#Grid refinement method and parameters
grid_refinement
// Available methods:
// double_points: []
// easy: [parameter]
// grcar_smooke: [param1, param2, param3] (legacy token: grcarsmooke)
// pearson: [param1, param2]
// twopnt: [param1, param2, param3]
pearson: [0.1, 1.5]

# Output and visualization options
postprocessing
# Generate gnuplot output
gnuplot: true
# Save results to CSV
save_to_csv: false
# Output filename prefix
filename: output_name
"#;

    let mut file = File::create("template.txt").map_err(|err| {
        crate::ReactorsBVP::SimpleReactorBVP::ReactorError::CalculationError(format!(
            "Failed to create template.txt: {}",
            err
        ))
    })?;
    file.write_all(template_content.as_bytes()).map_err(|err| {
        crate::ReactorsBVP::SimpleReactorBVP::ReactorError::CalculationError(format!(
            "Failed to write template.txt: {}",
            err
        ))
    })?;
    info!("Template created: template.txt");
    Ok(())
}
////////////////////////////////////////////////////TESTS///////////////////////////////////////////////////
#[cfg(test)]
mod tests {
    use super::*;
    use crate::ReactorsBVP::SimpleReactorBVP::ReactorError;
    use std::fs::File;
    use std::io::Write;
    use tempfile::tempdir;
    const task_content: &str = r#"
        initial_guess
        universal:1e-2
        process_conditions
        problem_name: Some(HMXTest)
        problem_description: Some(HMXdecompositiontest)
        substances: HMX, HMXprod
        t0: 0.0
        t_end: 1.0
        n_steps: 25
        arg:x
        Tm: 1500.0
        L: 5e-4
        dT: 600.0
        T_scale: 600.0
        P: 1e6
        Cp: 1464.4
        Lambda: 0.07
        m: 0.0043
        M: 0.0342
        thermal_effects: [102000.0]
        groups:true
        boundary_condition
        HMX: 0.999
        HMXprod: 0.001
        T: 800.0
        diffusion_coefficients
        HMX: 0.000009296
        HMXprod: 0.000009296
        HMX
        H: 4
        N: 8
        C: 8
        O: 8
        HMXprod
        H: 6
        C: 1
        O: 1
        reactions
        HMX=>10HMXprod: [130000.0, 0.0, 20920.0, 102000.0]
        solver_settings
        scheme: forward
        method: Sparse
        strategy: Damped
        linear_sys_method: None
        abs_tolerance: 1e-7
        max_iterations: 100
        loglevel: Some(info)
        dont_save_logs: true
        bounds
        C: -10.0, 10.0
        J:  -1e20, 1e20
        Teta:-100.0, 100.0
        q: -1e20, 1e20
        rel_tolerance
        C: 1e-7
        J: 1e-7
        Teta: 1e-7
        q:  1e-7
        strategy_params
        max_jac: Some(3)
        max_damp_iter: Some(10)
        damp_factor: Some(0.5)
        # Adaptive grid refinement settings (optional)
        adaptive_strategy
        # Refinement version
        version: 1
        # Maximum refinement iterations
        max_refinements: 3
                
        #Grid refinement method and parameters
        grid_refinement
        // Available methods:
        // double_points: []
        // easy: [parameter]
        // grcar_smooke: [param1, param2, param3] (legacy token: grcarsmooke)
        // pearson: [param1, param2]
        // twopnt: [param1, param2, param3]
        grcar_smooke: [0.05, 0.05, 1.25]
        postprocessing
        gnuplot:true
        save_to_csv:false
        filename: meow
        "#;
    #[test]
    fn parse_directly() {
        let mut parser = DocumentParser::new(task_content.to_string());
        let result: DocumentMap = parser.parse_document().unwrap().clone();
        info!("result {:?}", result);
    }
    #[test]
    fn check_substances() {
        let mut parser = DocumentParser::new(task_content.to_string());
        let result: DocumentMap = parser.parse_document().unwrap().clone();
        let process_conditions = result.get("process_conditions").unwrap();
        let substances = process_conditions
            .get("substances")
            .unwrap()
            .clone()
            .unwrap();

        info!("substances {:?}", substances[0]);
    }

    #[test]
    fn test_solve_from_parsed() {
        let mut parser = DocumentParser::new(task_content.to_string());
        let result: DocumentMap = parser.parse_document().unwrap().clone();
        let mut reactor = SimpleReactorTask::new();
        reactor
            .solve_from_parsed(result)
            .expect("Failed to solve parsed task");
        assert!(reactor.solver.solution.is_some());
        assert!(reactor.solver.x_mesh.is_some());
        assert!(!reactor.solver.x_mesh.as_ref().unwrap().is_empty());
        info!("solution {:?}", reactor.solver.solution);
    }
    #[test]
    fn parsng_task_elementary() {
        let mut reactor = SimpleReactorTask::new();

        // Create temporary directory
        let temp_dir = tempdir().expect("Failed to create temp dir");
        let file_path = temp_dir.path().join("hmx_task.txt");

        // HMX data from create_hmx function

        // Write to temporary file
        let mut file = File::create(&file_path).expect("Failed to create temp file");
        file.write_all(task_content.as_bytes())
            .expect("Failed to write to temp file");

        // Test parsing
        let mut parser = reactor
            .parse_file(Some(file_path))
            .expect("Failed to parse file");
        let _ = parser.parse_document();
        info!("parser {:?}", parser);

        reactor
            .parse_parameters_with_exact_names(&mut parser)
            .expect("Failed to parse parameters");

        // Verify parsed data
        assert_eq!(reactor.problem_name, Some("HMXTest".to_string()));
        assert_eq!(
            reactor.kindata.substances,
            vec!["HMX".to_string(), "HMXprod".to_string()]
        );

        assert_eq!(reactor.Tm, 1500.0);
        assert_eq!(reactor.P, 1e6);
        assert_eq!(reactor.Cp, 1464.4);
        assert_eq!(reactor.Lambda, 0.07);
        assert_eq!(reactor.m, 0.0043);
        assert_eq!(reactor.M, 0.0342);
        assert_eq!(reactor.thermal_effects, vec![102000.0]);

        // Verify boundary conditions
        assert_eq!(reactor.boundary_condition.get("HMX"), Some(&0.999));
        assert_eq!(reactor.boundary_condition.get("HMXprod"), Some(&0.001));
        assert_eq!(reactor.boundary_condition.get("T"), Some(&800.0));

        // Verify diffusion coefficients
        assert_eq!(reactor.Diffusion.get("HMX"), Some(&0.000009296));
        assert_eq!(reactor.Diffusion.get("HMXprod"), Some(&0.000009296));

        // Verify groups
        assert!(reactor.kindata.groups.is_some());
        let groups = reactor.kindata.groups.as_ref().unwrap();
        assert_eq!(groups.get("HMX").unwrap().get("H"), Some(&4));
        assert_eq!(groups.get("HMX").unwrap().get("N"), Some(&8));
        assert_eq!(groups.get("HMXprod").unwrap().get("H"), Some(&6));
        assert_eq!(groups.get("HMXprod").unwrap().get("C"), Some(&1));

        // Verify reactions were parsed
        assert_eq!(reactor.kindata.vec_of_equations.len(), 1);
        assert_eq!(reactor.kindata.vec_of_equations[0], "HMX=>10HMXprod");
    }

    #[test]
    fn test_parse_toleranse_and_bounds() {
        let mut reactor = SimpleReactorTask::new();

        // Create temporary directory
        let temp_dir = tempdir().expect("Failed to create temp dir");
        let file_path = temp_dir.path().join("test_bounds.txt");

        // Write to temporary file
        let mut file = File::create(&file_path).expect("Failed to create temp file");
        file.write_all(task_content.as_bytes())
            .expect("Failed to write to temp file");

        // Parse file
        let mut parser = reactor
            .parse_file(Some(file_path))
            .expect("Failed to parse file");
        let _ = parser.parse_document();

        // Test parsing tolerances and bounds
        let (tolerance_config, bounds_config) = reactor
            .parse_toleranse_and_bounds(&mut parser)
            .expect("Failed to parse tolerances and bounds");

        // Verify tolerance config
        assert_eq!(tolerance_config.C, 1e-7);
        assert_eq!(tolerance_config.J, 1e-7);
        assert_eq!(tolerance_config.Teta, 1e-7);
        assert_eq!(tolerance_config.q, 1e-7);

        // Verify bounds config
        assert_eq!(bounds_config.C, (-10.0, 10.0));
        assert_eq!(bounds_config.J, (-1e20, 1e20));
        assert_eq!(bounds_config.Teta, (-100.0, 100.0));
        assert_eq!(bounds_config.q, (-1e20, 1e20));
    }

    #[test]
    fn test_parse_solver_settings_accepts_aot_backend_options() {
        let mut solver_settings = HashMap::new();
        solver_settings.insert(
            "scheme".to_string(),
            Some(vec![Value::String("forward".to_string())]),
        );
        solver_settings.insert(
            "strategy".to_string(),
            Some(vec![Value::String("Damped".to_string())]),
        );
        solver_settings.insert(
            "method".to_string(),
            Some(vec![Value::String("Sparse".to_string())]),
        );
        solver_settings.insert("abs_tolerance".to_string(), Some(vec![Value::Float(1e-8)]));
        solver_settings.insert("max_iterations".to_string(), Some(vec![Value::Usize(50)]));
        solver_settings.insert(
            "generated_backend".to_string(),
            Some(vec![Value::String("banded_aot_tcc".to_string())]),
        );
        solver_settings.insert(
            "aot_compiler".to_string(),
            Some(vec![Value::String("zig".to_string())]),
        );
        solver_settings.insert(
            "aot_build_policy".to_string(),
            Some(vec![Value::String("require_prebuilt".to_string())]),
        );
        solver_settings.insert(
            "aot_compile_preset".to_string(),
            Some(vec![Value::String("fast_build".to_string())]),
        );
        solver_settings.insert(
            "aot_execution_policy".to_string(),
            Some(vec![Value::String("parallel".to_string())]),
        );
        solver_settings.insert("aot_target_chunks".to_string(), Some(vec![Value::Usize(4)]));

        let mut document = DocumentMap::new();
        document.insert("solver_settings".to_string(), solver_settings);

        let parsed =
            parse_solver_settings_from_map(&document).expect("AOT solver settings should parse");
        let options = parsed
            .solver_backend_config
            .to_rusted_options()
            .expect("Parsed AOT settings should build RustedSciThe options");

        assert_eq!(
            options.generated_backend_config.aot_codegen_backend,
            RustedSciThe::symbolic::codegen::codegen_aot_driver::AotCodegenBackend::Zig
        );
        assert_eq!(options.generated_backend_config.aot_c_compiler, None);
        assert!(matches!(
            options.generated_backend_config.aot_build_policy,
            RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::AotBuildPolicy::RequirePrebuilt
        ));
        assert!(matches!(
            options.generated_backend_config.aot_execution_policy,
            RustedSciThe::numerical::BVP_Damp::generated_solver_handoff::AotExecutionPolicy::Parallel(_)
        ));
        assert!(
            options
                .generated_backend_config
                .aot_chunking_policy
                .residual
                .is_some()
        );
        assert!(
            options
                .generated_backend_config
                .aot_chunking_policy
                .sparse_jacobian
                .is_some()
        );
    }

    #[test]
    fn test_parse_tolerance_and_bounds_alias_matches_legacy_spelling() {
        let mut reactor = SimpleReactorTask::new();

        let temp_dir = tempdir().expect("Failed to create temp dir");
        let file_path = temp_dir.path().join("test_bounds_alias.txt");

        let mut file = File::create(&file_path).expect("Failed to create temp file");
        file.write_all(task_content.as_bytes())
            .expect("Failed to write to temp file");

        let mut legacy_parser = reactor
            .parse_file(Some(file_path.clone()))
            .expect("Failed to parse file");
        let _ = legacy_parser.parse_document();

        let mut alias_parser = reactor
            .parse_file(Some(file_path))
            .expect("Failed to parse file");
        let _ = alias_parser.parse_document();

        let legacy = reactor
            .parse_toleranse_and_bounds(&mut legacy_parser)
            .expect("Legacy spelling should still work");
        let alias = reactor
            .parse_tolerance_and_bounds(&mut alias_parser)
            .expect("New spelling should work too");

        assert_eq!(legacy.0.C, alias.0.C);
        assert_eq!(legacy.0.J, alias.0.J);
        assert_eq!(legacy.0.Teta, alias.0.Teta);
        assert_eq!(legacy.0.q, alias.0.q);
        assert_eq!(legacy.1.C, alias.1.C);
        assert_eq!(legacy.1.J, alias.1.J);
        assert_eq!(legacy.1.Teta, alias.1.Teta);
        assert_eq!(legacy.1.q, alias.1.q);
    }

    #[test]
    fn test_parse_grid_refinement_alias_matches_legacy_spelling() {
        let mut reactor = SimpleReactorTask::new();

        let legacy_task = task_content.to_string();
        let alias_task = task_content
            .replace("grcarsmooke", "grcar_smooke")
            .replace("doubleoints", "double_points");

        let mut legacy_parser = DocumentParser::new(legacy_task);
        legacy_parser
            .parse_document()
            .expect("Legacy spelling should parse");
        let mut alias_parser = DocumentParser::new(alias_task);
        alias_parser
            .parse_document()
            .expect("Corrected spelling should parse");

        let legacy = reactor.parse_parameters_with_exact_names(&mut legacy_parser);
        let alias = reactor.parse_parameters_with_exact_names(&mut alias_parser);

        assert!(legacy.is_ok());
        assert!(alias.is_ok());
    }

    #[test]
    fn test_parse_toleranse_and_bounds_accepts_scalar_pairs() {
        let mut bounds = HashMap::new();
        bounds.insert(
            "C".to_string(),
            Some(vec![Value::Float(-10.0), Value::Float(10.0)]),
        );
        bounds.insert(
            "J".to_string(),
            Some(vec![Value::Float(-1e20), Value::Float(1e20)]),
        );
        bounds.insert(
            "Teta".to_string(),
            Some(vec![Value::Float(-100.0), Value::Float(100.0)]),
        );
        bounds.insert(
            "q".to_string(),
            Some(vec![Value::Float(-1e20), Value::Float(1e20)]),
        );

        let mut rel_tolerance = HashMap::new();
        rel_tolerance.insert("C".to_string(), Some(vec![Value::Float(1e-7)]));
        rel_tolerance.insert("J".to_string(), Some(vec![Value::Float(1e-7)]));
        rel_tolerance.insert("Teta".to_string(), Some(vec![Value::Float(1e-7)]));
        rel_tolerance.insert("q".to_string(), Some(vec![Value::Float(1e-7)]));

        let mut document = DocumentMap::new();
        document.insert("bounds".to_string(), bounds);
        document.insert("rel_tolerance".to_string(), rel_tolerance);

        let (tolerance_config, bounds_config) = parse_tolerance_and_bounds_from_map(&document)
            .expect("Failed to parse scalar-pair bounds");

        assert_eq!(tolerance_config.C, 1e-7);
        assert_eq!(bounds_config.C, (-10.0, 10.0));
        assert_eq!(bounds_config.J, (-1e20, 1e20));
        assert_eq!(bounds_config.Teta, (-100.0, 100.0));
        assert_eq!(bounds_config.q, (-1e20, 1e20));
    }

    #[test]
    fn test_validate_solver_output_matrix_rejects_non_finite_values() {
        let matrix = DMatrix::from_vec(2, 2, vec![1.0, 2.0, 3.0, f64::NAN]);
        let result = validate_solver_output_matrix(&matrix);
        assert!(result.is_err());
        match result {
            Err(ReactorError::InvalidNumericValue(msg)) => {
                assert!(msg.contains("non-finite value"));
            }
            _ => panic!("Expected InvalidNumericValue error"),
        }
    }

    #[test]
    fn test_parse_parameters_rejects_non_finite_pressure() {
        let mut reactor = SimpleReactorTask::new();
        let task = task_content.replace("P: 1e6", "P: NaN");
        let mut parser = DocumentParser::new(task);
        let result = reactor.parse_parameters_with_exact_names(&mut parser);
        assert!(matches!(result, Err(ReactorError::InvalidNumericValue(_))));
    }

    #[test]
    fn test_parse_parameters_with_exact_names_is_idempotent_for_preparsed_document() {
        let mut reactor = SimpleReactorTask::new();
        let mut parser = DocumentParser::new(task_content.to_string());
        parser.parse_document().expect("Failed to parse document");

        reactor
            .parse_parameters_with_exact_names(&mut parser)
            .expect("First parse should succeed");
        reactor
            .parse_parameters_with_exact_names(&mut parser)
            .expect("Second parse should also succeed");

        assert_eq!(reactor.kindata.substances, vec!["HMX", "HMXprod"]);
        assert_eq!(reactor.kindata.vec_of_equations, vec!["HMX=>10HMXprod"]);
        assert_eq!(reactor.boundary_condition.get("HMX"), Some(&0.999));
        assert_eq!(reactor.boundary_condition.get("HMXprod"), Some(&0.001));
    }

    #[test]
    fn test_setup_bvp_from_file() {
        let mut reactor = SimpleReactorTask::new();
        // Create temporary directory
        let temp_dir = tempdir().expect("Failed to create temp dir");
        let file_path = temp_dir.path().join("hmx_task.txt");
        // Write to temporary file
        let mut file = File::create(&file_path).expect("Failed to create temp file");
        file.write_all(task_content.as_bytes())
            .expect("Failed to write to temp file");

        // Test parsing
        let mut parser = reactor
            .parse_file(Some(file_path))
            .expect("Failed to parse file");
        let _ = parser.parse_document();
        // println!("parser {:?}", parser);

        reactor
            .parse_parameters_with_exact_names(&mut parser)
            .expect("Failed to parse parameters");

        reactor.setup_bvp().expect("Failed to setup BVP");
        let rates = reactor.map_eq_rate.clone();
        for (eq, rate) in rates {
            info!("reaction {} rate {}", eq, rate);
        }
        info!("\n \n");
        let system = reactor.map_of_equations.clone();
        for (subs, (variable, eq)) in system {
            info!("subs: {} | variable: {} | eq: {} | \n", subs, variable, eq);
        }
        let bc = &reactor.solver.BorderConditions;
        info!("bc {:?}", bc);
        info!(" unknowns{:?}", reactor.solver.unknowns);
    }

    #[test]
    fn test_solve_from_map() {
        let mut reactor = SimpleReactorTask::new();
        let temp_dir = tempdir().expect("Failed to create temp dir");
        let file_path = temp_dir.path().join("hmx_task.txt");
        // Write to temporary file
        let mut file = File::create(&file_path).expect("Failed to create temp file");
        file.write_all(task_content.as_bytes())
            .expect("Failed to write to temp file");

        // Test parsing
        let mut parser = reactor
            .parse_file(Some(file_path))
            .expect("Failed to parse file");
        let _ = parser.parse_document();
        reactor
            .solve_from_map(parser)
            .expect("Failed to solve from map");
        // println!("parser {:?}", parser);
    }
}