use rand::Rng;
use rand_distr::{Normal, Distribution};
use ndarray::Array1;

/// Computes the thermal conductivity for a given temperature using a linear model.
///
/// This function models the thermal conductivity as a linear function of temperature,
/// where the conductivity is calculated as the product of a coefficient and the temperature,
/// plus an intercept. This simple model serves as a basis for demonstrating uncertainty propagation.
///
/// # Arguments
///
/// * `temp` - A f64 value representing the temperature.
/// * `coeff` - A f64 value representing the temperature coefficient in the model.
/// * `intercept` - A f64 value representing the baseline thermal conductivity.
///
/// # Returns
///
/// A f64 value representing the predicted thermal conductivity.
fn thermal_conductivity_model(temp: f64, coeff: f64, intercept: f64) -> f64 {
    coeff * temp + intercept
}

/// Simulates thermal conductivity data by incorporating uncertainty in temperature measurements.
///
/// For each temperature in the input slice, this function adds Gaussian noise to simulate measurement uncertainty
/// and then computes the corresponding thermal conductivity using the linear model.
///
/// # Arguments
///
/// * `temps` - A slice of f64 values representing measured temperatures.
/// * `true_coeff` - A f64 value representing the true coefficient for temperature dependence.
/// * `true_intercept` - A f64 value representing the true baseline thermal conductivity.
/// * `temp_noise` - A f64 value representing the standard deviation of the noise added to the temperature.
///
/// # Returns
///
/// A vector of f64 values representing the simulated thermal conductivity data.
fn simulate_conductivity_data(
    temps: &[f64],
    true_coeff: f64,
    true_intercept: f64,
    temp_noise: f64
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    // Define a normal distribution for the temperature noise.
    let noise_dist = Normal::new(0.0, temp_noise).unwrap();
    // For each temperature, add noise and compute the conductivity.
    temps.iter().map(|&t| {
        let noisy_temp = t + noise_dist.sample(&mut rng);
        thermal_conductivity_model(noisy_temp, true_coeff, true_intercept)
    }).collect()
}

/// Quantifies the uncertainty in conductivity predictions using Monte Carlo sampling.
///
/// This function estimates the uncertainty in the thermal conductivity predictions by sampling from the distributions
/// of both the model parameters (coefficient and intercept) and the temperature measurements. For each temperature,
/// it randomly selects a parameter sample and adds noise to the temperature before computing the conductivity.
/// The result is a vector of predictions that reflects the propagated uncertainty from all input sources.
///
/// # Arguments
///
/// * `temps` - A slice of f64 values representing the temperature data.
/// * `coeff_samples` - A slice of f64 values representing sampled values for the coefficient parameter.
/// * `intercept_samples` - A slice of f64 values representing sampled values for the intercept parameter.
/// * `temp_noise` - A f64 value representing the standard deviation of the temperature measurement noise.
///
/// # Returns
///
/// A vector of f64 values representing the predicted thermal conductivities with uncertainty.
fn quantify_uncertainty(
    temps: &[f64],
    coeff_samples: &[f64],
    intercept_samples: &[f64],
    temp_noise: f64
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    let noise_dist = Normal::new(0.0, temp_noise).unwrap();

    // Initialize a vector to hold the predictions.
    let mut predictions = vec![0.0; temps.len()];
    for (&t, prediction) in temps.iter().zip(predictions.iter_mut()) {
        // For each temperature, sample a noisy temperature.
        let temp_sample = t + noise_dist.sample(&mut rng);
        // Randomly select a coefficient and intercept sample.
        let coeff_sample = coeff_samples[rng.gen_range(0..coeff_samples.len())];
        let intercept_sample = intercept_samples[rng.gen_range(0..intercept_samples.len())];
        // Compute the predicted conductivity using the sampled parameters and temperature.
        *prediction = thermal_conductivity_model(temp_sample, coeff_sample, intercept_sample);
    }
    predictions
}

/// Main function demonstrating uncertainty quantification in a material property prediction model.
///
/// In this example, we simulate a set of temperature measurements, generate corresponding thermal conductivity data
/// using a linear model with known true parameters, and then perform uncertainty quantification using Monte Carlo
/// sampling. Parameter samples for the model are generated around the true values with small random variations to
/// simulate uncertainty in parameter estimation. The final predictions, which include propagated uncertainties, are
/// printed alongside their corresponding temperatures.
fn main() {
    // Define the temperature data as measured input values (in Kelvin).
    let temps = vec![300.0, 350.0, 400.0, 450.0, 500.0];

    // Define the true parameters for the thermal conductivity model.
    let true_coeff = 0.02;    // True coefficient for temperature dependence.
    let true_intercept = 1.0;   // True baseline thermal conductivity.
    let temp_noise = 5.0;       // Standard deviation of the noise in temperature measurements.

    // Simulate the thermal conductivity data by incorporating measurement noise.
    let data = simulate_conductivity_data(&temps, true_coeff, true_intercept, temp_noise);

    // Generate samples for the model parameters to simulate uncertainty in parameter estimation.
    let coeff_samples: Vec<f64> = (0..100)
        .map(|_| true_coeff + rand::thread_rng().gen_range(-0.001..0.001))
        .collect();
    let intercept_samples: Vec<f64> = (0..100)
        .map(|_| true_intercept + rand::thread_rng().gen_range(-0.05..0.05))
        .collect();

    // Use Monte Carlo sampling to propagate uncertainty in both the temperature data and the model parameters.
    let predictions = quantify_uncertainty(&temps, &coeff_samples, &intercept_samples, temp_noise);

    // Display the predictions along with the corresponding temperature values.
    for (temp, prediction) in temps.iter().zip(predictions.iter()) {
        println!("Temperature: {:.1} K, Predicted Conductivity: {:.4} W/mK", temp, prediction);
    }
}

use rand::Rng;
use rand_distr::{Normal, Distribution};
use ndarray::Array1;

/// Defines a simplified fluid dynamics model where the fluid velocity is modeled as the ratio
/// of a pressure drop to the fluid's viscosity. This model assumes a direct inverse relationship
/// between viscosity and velocity.
///
/// # Arguments
///
/// * `pressure_drop` - A f64 value representing the pressure difference driving the fluid flow.
/// * `viscosity` - A f64 value representing the viscosity of the fluid.
///
/// # Returns
///
/// A f64 value representing the predicted fluid velocity.
fn fluid_velocity_model(pressure_drop: f64, viscosity: f64) -> f64 {
    pressure_drop / viscosity
}

/// Simulates fluid velocity data using Monte Carlo sampling.
///
/// This function generates multiple samples of fluid velocity by drawing random samples from a normal
/// distribution representing the uncertain viscosity. The pressure drop is treated as a fixed input.
/// The resulting collection of fluid velocities reflects the propagated uncertainty from the viscosity.
///
/// # Arguments
///
/// * `pressure_drop` - A f64 value representing the fixed pressure drop in the system.
/// * `viscosity_mean` - A f64 value representing the mean of the viscosity distribution.
/// * `viscosity_std` - A f64 value representing the standard deviation of the viscosity distribution.
/// * `num_samples` - The number of Monte Carlo samples to generate.
///
/// # Returns
///
/// A vector of f64 values containing the simulated fluid velocities.
fn monte_carlo_simulation(
    pressure_drop: f64,
    viscosity_mean: f64,
    viscosity_std: f64,
    num_samples: usize,
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    let viscosity_dist = Normal::new(viscosity_mean, viscosity_std).unwrap();

    // Generate samples by drawing from the viscosity distribution and computing the corresponding velocity.
    (0..num_samples)
        .map(|_| {
            let viscosity_sample = viscosity_dist.sample(&mut rng);
            fluid_velocity_model(pressure_drop, viscosity_sample)
        })
        .collect()
}

/// Computes basic statistics (mean and standard deviation) from a set of Monte Carlo samples.
///
/// This function aggregates the results of the simulation to provide an estimate of the expected fluid velocity
/// along with a measure of its uncertainty.
///
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the simulated fluid velocities.
///
/// # Returns
///
/// A tuple (mean, standard deviation) representing the central tendency and spread of the fluid velocity.
fn compute_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / (samples.len() as f64);
    (mean, variance.sqrt())
}

fn main() {
    // Define simulation parameters: pressure drop, and viscosity modeled as a normal distribution.
    let pressure_drop = 10.0; // Pressure drop in Pascals.
    let viscosity_mean = 1.0; // Mean viscosity in Pa·s.
    let viscosity_std = 0.2;  // Standard deviation of viscosity.
    let num_samples = 1000;   // Number of Monte Carlo samples.

    // Run Monte Carlo simulation to generate fluid velocity samples.
    let velocity_samples = monte_carlo_simulation(pressure_drop, viscosity_mean, viscosity_std, num_samples);

    // Compute the mean fluid velocity and its uncertainty from the simulation.
    let (mean_velocity, std_velocity) = compute_statistics(&velocity_samples);

    // Output the computed statistics.
    println!("Mean fluid velocity: {:.2} m/s", mean_velocity);
    println!("Uncertainty (standard deviation): {:.2} m/s", std_velocity);
}

use rand_distr::{Normal, Distribution};
use rand::thread_rng;
use rand::Rng; // Added import to bring the `gen` method into scope.
use ndarray::Array1;

/// Computes the probability density function (pdf) of a Gaussian distribution with mean `mu` and standard deviation `sigma` at the given value `x`.
///
/// This helper function implements the standard Gaussian pdf formula, which is used to evaluate the probability of a residual under a Gaussian noise model.
///
/// # Arguments
///
/// * `x` - The value at which the density is evaluated.
/// * `mu` - The mean of the distribution.
/// * `sigma` - The standard deviation of the distribution.
///
/// # Returns
///
/// A f64 value representing the probability density at `x`.
fn normal_pdf(x: f64, mu: f64, sigma: f64) -> f64 {
    let coeff = 1.0 / ((2.0 * std::f64::consts::PI).sqrt() * sigma);
    let exponent = -((x - mu).powi(2)) / (2.0 * sigma * sigma);
    coeff * exponent.exp()
}

/// Defines a simplified fluid dynamics model where the fluid velocity is modeled as the ratio
/// of a pressure drop to the fluid's viscosity. This model assumes a direct inverse relationship
/// between viscosity and velocity.
///
/// # Arguments
///
/// * `pressure_drop` - A f64 value representing the pressure difference driving the fluid flow.
/// * `viscosity` - A f64 value representing the viscosity of the fluid.
///
/// # Returns
///
/// A f64 value representing the predicted fluid velocity.
fn fluid_velocity_model(pressure_drop: f64, viscosity: f64) -> f64 {
    pressure_drop / viscosity
}

/// Simulates fluid velocity data using Monte Carlo sampling.
///
/// This function generates multiple samples of fluid velocity by drawing random samples from a normal
/// distribution representing the uncertain viscosity. The pressure drop is treated as a fixed input,
/// and the resulting collection of fluid velocities reflects the propagated uncertainty from the viscosity.
///
/// # Arguments
///
/// * `pressure_drop` - A f64 value representing the fixed pressure drop in the system.
/// * `viscosity_mean` - A f64 value representing the mean of the viscosity distribution.
/// * `viscosity_std` - A f64 value representing the standard deviation of the viscosity distribution.
/// * `num_samples` - The number of Monte Carlo samples to generate.
///
/// # Returns
///
/// A vector of f64 values containing the simulated fluid velocities.
fn monte_carlo_simulation(
    pressure_drop: f64,
    viscosity_mean: f64,
    viscosity_std: f64,
    num_samples: usize,
) -> Vec<f64> {
    let mut rng = thread_rng();
    let viscosity_dist = Normal::new(viscosity_mean, viscosity_std).unwrap();

    // Generate samples by drawing from the viscosity distribution and computing the corresponding velocity.
    (0..num_samples)
        .map(|_| {
            let viscosity_sample = viscosity_dist.sample(&mut rng);
            fluid_velocity_model(pressure_drop, viscosity_sample)
        })
        .collect()
}

/// Computes basic statistics (mean and standard deviation) from a set of Monte Carlo samples.
///
/// This function aggregates the simulation results to provide an estimate of the expected fluid velocity
/// along with a measure of its uncertainty.
///
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the simulated fluid velocities.
///
/// # Returns
///
/// A tuple (mean, standard deviation) representing the central tendency and spread of the fluid velocity.
fn compute_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / (samples.len() as f64);
    (mean, variance.sqrt())
}

/// Likelihood function that computes the probability of observing the data given a certain viscosity.
///
/// The likelihood is defined under a Gaussian noise model, where the difference between the observed fluid velocity
/// and the predicted velocity (computed by the fluid dynamics model) is assumed to be normally distributed.
///
/// # Arguments
///
/// * `data` - A slice of f64 values representing observed fluid velocities.
/// * `pressure_drop` - A f64 value representing the fixed pressure drop in the system.
/// * `viscosity` - A f64 value representing the viscosity parameter to evaluate.
/// * `noise` - A f64 value representing the standard deviation of the measurement noise.
///
/// # Returns
///
/// A f64 value representing the product of the probability densities for all observed data points.
fn likelihood(data: &[f64], pressure_drop: f64, viscosity: f64, noise: f64) -> f64 {
    data.iter().map(|&v| {
        let predicted_velocity = fluid_velocity_model(pressure_drop, viscosity);
        normal_pdf(v - predicted_velocity, 0.0, noise)
    }).product()
}

/// Performs a Bayesian update to estimate the viscosity parameter using the Metropolis-Hastings algorithm.
///
/// This function iteratively proposes new values for the viscosity parameter and accepts or rejects these proposals
/// based on the likelihood of the observed data under the current and proposed values. The process yields an updated
/// estimate that reflects both the prior belief and the observed data.
///
/// # Arguments
///
/// * `data` - A slice of f64 values representing observed fluid velocities.
/// * `pressure_drop` - A f64 value representing the fixed pressure drop.
/// * `prior_mean` - A f64 value representing the prior mean for the viscosity.
/// * `prior_std` - A f64 value representing the prior standard deviation for the viscosity (not used explicitly here).
/// * `noise` - A f64 value representing the standard deviation of the measurement noise.
///
/// # Returns
///
/// A f64 value representing the updated viscosity estimate after performing the Bayesian update.
fn bayesian_update(
    data: &[f64],
    pressure_drop: f64,
    prior_mean: f64,
    _prior_std: f64,
    noise: f64,
) -> f64 {
    let mut rng = thread_rng();
    let proposal_dist = Normal::new(0.0, 0.1).unwrap();

    // Initialize the viscosity parameter with the prior mean.
    let mut current_viscosity = prior_mean;

    // Run the Metropolis-Hastings algorithm for a fixed number of iterations.
    for _ in 0..1000 {
        let proposed_viscosity = current_viscosity + proposal_dist.sample(&mut rng);
        let likelihood_current = likelihood(data, pressure_drop, current_viscosity, noise);
        let likelihood_proposed = likelihood(data, pressure_drop, proposed_viscosity, noise);

        // Compute the acceptance ratio.
        let acceptance_ratio = likelihood_proposed / likelihood_current;

        // Accept the proposed viscosity with probability equal to the acceptance ratio.
        if rng.gen::<f64>() < acceptance_ratio {
            current_viscosity = proposed_viscosity;
        }
    }

    current_viscosity
}

/// Main function demonstrating a probabilistic approach to uncertainty quantification using Bayesian inference.
///
/// In this example, observed fluid velocity data is simulated using a Monte Carlo method applied to a simple fluid
/// dynamics model with uncertain viscosity. The Metropolis-Hastings algorithm is then used to update our estimate
/// of the viscosity parameter, incorporating the uncertainty from the observed data. The updated viscosity estimate,
/// along with basic statistics of the simulated velocities, provides insight into the uncertainty propagation in the model.
fn main() {
    // Define simulation parameters.
    let pressure_drop = 10.0; // Pressure drop in Pascals.
    let true_viscosity = 1.0; // True viscosity in Pa·s.
    let noise = 0.1;        // Measurement noise standard deviation.
    let num_samples_mc = 10; // A small number of samples for demonstration purposes.

    // Simulate observed fluid velocity data using Monte Carlo sampling.
    let data = monte_carlo_simulation(pressure_drop, true_viscosity, noise, num_samples_mc);

    // Define prior beliefs for viscosity.
    let prior_mean = 1.0;
    let prior_std = 0.5; // Although not used directly, this can inform more advanced priors.

    // Update the viscosity estimate using Bayesian inference.
    let estimated_viscosity = bayesian_update(&data, pressure_drop, prior_mean, prior_std, noise);

    // Output the estimated viscosity.
    println!("Estimated viscosity: {:.2} Pa·s", estimated_viscosity);

    // Compute and display basic statistics for the Monte Carlo simulation.
    let (mean_velocity, std_velocity) = compute_statistics(&data);
    println!("Monte Carlo simulated velocity: mean = {:.2}, std = {:.2}", mean_velocity, std_velocity);
}

use ndarray::Array1;

/// Computes the structural stress given a load and a cross-sectional area.
/// This simplified model assumes that stress is defined as the ratio of load to area.
///
/// # Arguments
///
/// * `load` - A f64 value representing the applied load (in Newtons).
/// * `area` - A f64 value representing the cross-sectional area (in square meters).
///
/// # Returns
///
/// A f64 value representing the computed stress (in Pascals).
fn structural_stress(load: f64, area: f64) -> f64 {
    load / area
}

/// Performs interval analysis for uncertain load and area values.
/// The analysis computes the minimum and maximum possible stress by evaluating the model at the endpoints
/// of the input intervals. The minimum stress is obtained by taking the lowest load and the highest area,
/// while the maximum stress is obtained by taking the highest load and the lowest area.
///
/// # Arguments
///
/// * `load_interval` - A tuple (f64, f64) representing the lower and upper bounds of the load.
/// * `area_interval` - A tuple (f64, f64) representing the lower and upper bounds of the cross-sectional area.
///
/// # Returns
///
/// A tuple (min_stress, max_stress) representing the range of possible stresses.
fn interval_analysis(load_interval: (f64, f64), area_interval: (f64, f64)) -> (f64, f64) {
    // Compute minimum stress using minimum load and maximum area.
    let min_stress = structural_stress(load_interval.0, area_interval.1);
    // Compute maximum stress using maximum load and minimum area.
    let max_stress = structural_stress(load_interval.1, area_interval.0);
    (min_stress, max_stress)
}

fn main() {
    // Define the uncertain input intervals for load and area.
    let load_interval = (1000.0, 1200.0); // Load in Newtons.
    let area_interval = (0.05, 0.1);        // Cross-sectional area in square meters.

    // Perform interval analysis to compute the range of possible stresses.
    let (min_stress, max_stress) = interval_analysis(load_interval, area_interval);

    // Output the computed stress range.
    println!("Stress range: {:.2} - {:.2} Pascals", min_stress, max_stress);
}

/// Computes the structural stress based on load and cross-sectional area.
/// This function assumes that stress is the ratio of load to area.
///
/// # Arguments
///
/// * `load` - A f64 value representing the load in Newtons.
/// * `area` - A f64 value representing the cross-sectional area in square meters.
///
/// # Returns
///
/// A f64 value representing the calculated stress in Pascals.
fn structural_stress(load: f64, area: f64) -> f64 {
    load / area
}

/// Defines fuzzy membership functions for the load variable.
/// The function returns a tuple representing the degree of membership for the load being
/// classified as low, medium, and high respectively.
///
/// # Arguments
///
/// * `load` - A f64 value representing the load.
///
/// # Returns
///
/// A tuple (low, medium, high) where each value is between 0.0 and 1.0.
fn fuzzy_load(load: f64) -> (f64, f64, f64) {
    // Define membership for "low" load: full membership for loads <= 500,
    // decreasing linearly to 0 for loads between 500 and 1000.
    let low = if load <= 500.0 {
        1.0
    } else if load < 1000.0 {
        1.0 - (load - 500.0) / 500.0
    } else {
        0.0
    };

    // Define membership for "medium" load: increases from 0 at 500, reaches 1 at 1000,
    // and then decreases back to 0 at 1500.
    let medium = if load >= 500.0 && load <= 1500.0 {
        if load <= 1000.0 {
            (load - 500.0) / 500.0
        } else {
            1.0 - (load - 1000.0) / 500.0
        }
    } else {
        0.0
    };

    // Define membership for "high" load: 0 for loads < 1000, increasing linearly thereafter.
    let high = if load >= 1000.0 {
        (load - 1000.0) / 500.0
    } else {
        0.0
    };

    (low, medium, high)
}

/// Uses fuzzy logic to assess structural safety by estimating the stress based on imprecise load data.
/// This function applies fuzzy membership functions to the load, computes the corresponding stress,
/// and then defuzzifies the result using a weighted average of the stress values for each category.
///
/// # Arguments
///
/// * `load` - A f64 value representing the load in Newtons.
/// * `area` - A f64 value representing the cross-sectional area in square meters.
///
/// # Returns
///
/// A f64 value representing the estimated structural stress based on fuzzy logic.
fn fuzzy_inference(load: f64, area: f64) -> f64 {
    // Obtain fuzzy membership values for the load.
    let (low_load, medium_load, high_load) = fuzzy_load(load);

    // Compute stress using the structural model.
    let stress = structural_stress(load, area);

    // In a more complex model, different rules would apply; here we simply weight the computed stress
    // by the degrees of membership in each category.
    let weighted_stress = stress * (low_load + medium_load + high_load);
    let total_membership = low_load + medium_load + high_load;

    // Avoid division by zero in case of no membership (should not happen with properly defined fuzzy sets).
    if total_membership == 0.0 {
        stress
    } else {
        weighted_stress / total_membership
    }
}

fn main() {
    // Define a sample load and cross-sectional area.
    let load = 800.0; // Load in Newtons.
    let area = 0.08;  // Cross-sectional area in square meters.

    // Perform fuzzy inference to estimate the structural stress.
    let safety_stress = fuzzy_inference(load, area);

    // Output the estimated stress, providing insight into structural safety under imprecise measurements.
    println!("Estimated stress based on fuzzy logic: {:.2} Pascals", safety_stress);
}

use ndarray::Array1;

/// A simple weather prediction model that estimates temperature based on humidity and pressure.
/// In this simplified linear model, temperature is assumed to be a weighted sum of humidity and pressure.
///
/// # Arguments
///
/// * `humidity` - A f64 value representing humidity (in percentage).
/// * `pressure` - A f64 value representing atmospheric pressure (in hPa).
///
/// # Returns
///
/// A f64 value representing the predicted temperature.
fn weather_model(humidity: f64, pressure: f64) -> f64 {
    0.6 * humidity + 0.4 * pressure
}

/// Computes the local sensitivity of the weather model output with respect to the input parameters.
/// This is achieved by slightly perturbing each input and evaluating the change in the output,
/// thereby approximating the partial derivatives.
///
/// # Arguments
///
/// * `humidity` - The nominal value of humidity.
/// * `pressure` - The nominal value of pressure.
/// * `delta` - A small perturbation value to compute finite differences.
///
/// # Returns
///
/// A tuple (humidity_sensitivity, pressure_sensitivity) representing the approximate derivatives.
fn local_sensitivity(humidity: f64, pressure: f64, delta: f64) -> (f64, f64) {
    let base_output = weather_model(humidity, pressure);

    // Approximate the sensitivity with respect to humidity by increasing humidity by delta.
    let humidity_sensitivity = (weather_model(humidity + delta, pressure) - base_output) / delta;

    // Approximate the sensitivity with respect to pressure by increasing pressure by delta.
    let pressure_sensitivity = (weather_model(humidity, pressure + delta) - base_output) / delta;

    (humidity_sensitivity, pressure_sensitivity)
}

fn main() {
    // Define nominal values for humidity and pressure.
    let humidity = 70.0;  // Example humidity in percentage.
    let pressure = 1000.0; // Example pressure in hPa.

    // Set a small perturbation for the sensitivity calculation.
    let delta = 1.0;

    // Compute the local sensitivities of the weather model output.
    let (humidity_sensitivity, pressure_sensitivity) = local_sensitivity(humidity, pressure, delta);

    println!("Local sensitivity with respect to humidity: {:.4}", humidity_sensitivity);
    println!("Local sensitivity with respect to pressure: {:.4}", pressure_sensitivity);
}

use rand::Rng;
use rand_distr::{Normal, Distribution};

/// The same weather model as before, predicting temperature based on humidity and pressure.
///
/// # Arguments
///
/// * `humidity` - A f64 value representing humidity.
/// * `pressure` - A f64 value representing pressure.
///
/// # Returns
///
/// A f64 value representing the predicted temperature.
fn weather_model(humidity: f64, pressure: f64) -> f64 {
    0.6 * humidity + 0.4 * pressure
}

/// Generates random samples for humidity and pressure for global sensitivity analysis.
/// The samples are drawn from normal distributions centered at the given means with a specified standard deviation.
///
/// # Arguments
///
/// * `num_samples` - The number of samples to generate.
/// * `mean_humidity` - The mean value for humidity.
/// * `mean_pressure` - The mean value for pressure.
/// * `std_dev` - The standard deviation for both humidity and pressure.
///
/// # Returns
///
/// A vector of tuples, where each tuple contains a pair of sampled (humidity, pressure) values.
fn generate_samples(num_samples: usize, mean_humidity: f64, mean_pressure: f64, std_dev: f64) -> Vec<(f64, f64)> {
    let mut rng = rand::thread_rng();
    let humidity_dist = Normal::new(mean_humidity, std_dev).unwrap();
    let pressure_dist = Normal::new(mean_pressure, std_dev).unwrap();

    (0..num_samples)
        .map(|_| {
            let humidity_sample = humidity_dist.sample(&mut rng);
            let pressure_sample = pressure_dist.sample(&mut rng);
            (humidity_sample, pressure_sample)
        })
        .collect()
}

/// Computes Sobol indices to quantify the contribution of each input to the output variance.
/// This simplified implementation calculates the variance of the model output and the first-order effects
/// of each parameter by fixing the other parameter at its mean value.
///
/// # Arguments
///
/// * `samples` - A slice of tuples containing the sampled input parameters.
/// * `mean_output` - The mean output of the model computed from the samples.
///
/// # Returns
///
/// A tuple (sobol_humidity, sobol_pressure) representing the normalized contribution of humidity and pressure,
/// respectively, to the total output variance.
fn sobol_indices(samples: &[(f64, f64)], mean_output: f64) -> (f64, f64) {
    let num_samples = samples.len() as f64;
    
    // Compute the total variance of the model output.
    let total_variance: f64 = samples.iter()
        .map(|(humidity, pressure)| {
            let output = weather_model(*humidity, *pressure);
            (output - mean_output).powi(2)
        })
        .sum::<f64>() / num_samples;

    // Compute the first-order effect of humidity by varying humidity while fixing pressure at its mean.
    let humidity_variance: f64 = samples.iter()
        .map(|(humidity, _)| {
            // For this demonstration, assume the mean pressure is derived from the overall mean output.
            // In practice, the mean value of pressure should be used.
            let fixed_pressure = mean_output / 0.4;
            let output = weather_model(*humidity, fixed_pressure);
            (output - mean_output).powi(2)
        })
        .sum::<f64>() / num_samples;

    // Compute the first-order effect of pressure by varying pressure while fixing humidity at its mean.
    let pressure_variance: f64 = samples.iter()
        .map(|(_, pressure)| {
            // Similarly, assume the mean humidity is derived from the overall mean output.
            let fixed_humidity = mean_output / 0.6;
            let output = weather_model(fixed_humidity, *pressure);
            (output - mean_output).powi(2)
        })
        .sum::<f64>() / num_samples;

    // Normalize the variances to obtain Sobol indices.
    (humidity_variance / total_variance, pressure_variance / total_variance)
}

fn main() {
    // Set parameters for generating input samples.
    let num_samples = 1000;
    let mean_humidity = 70.0;   // Mean humidity in percentage.
    let mean_pressure = 1000.0; // Mean pressure in hPa.
    let std_dev = 10.0;         // Standard deviation for both inputs.

    // Generate random samples for humidity and pressure.
    let samples = generate_samples(num_samples, mean_humidity, mean_pressure, std_dev);

    // Compute the mean output of the weather model over the generated samples.
    let mean_output: f64 = samples.iter()
        .map(|(humidity, pressure)| weather_model(*humidity, *pressure))
        .sum::<f64>() / num_samples as f64;

    // Calculate Sobol indices to quantify the sensitivity of the model output to each input.
    let (humidity_sobol, pressure_sobol) = sobol_indices(&samples, mean_output);

    println!("Sobol index for humidity: {:.4}", humidity_sobol);
    println!("Sobol index for pressure: {:.4}", pressure_sobol);
}

use rand::Rng;
use rand::thread_rng; // May emit a deprecation warning; consider updating the API if desired.
use ndarray::{Array1, Array2, ArrayViewMut2};
use ndarray::array;

/// Inverts a square matrix using Gauss-Jordan elimination.
/// Returns Some(inverse) if successful, or None if the matrix is singular.
///
/// # Arguments
///
/// * `matrix` - An Array2<f64> representing the square matrix.
fn invert_matrix(matrix: &Array2<f64>) -> Option<Array2<f64>> {
    let n = matrix.nrows();
    if matrix.ncols() != n {
        return None; // Not a square matrix.
    }

    // Create an augmented matrix [A | I].
    let mut aug = Array2::<f64>::zeros((n, 2 * n));
    // Copy the original matrix A.
    for i in 0..n {
        for j in 0..n {
            aug[[i, j]] = matrix[[i, j]];
        }
    }
    // Set the right half to the identity matrix.
    for i in 0..n {
        aug[[i, n + i]] = 1.0;
    }

    // Perform Gauss-Jordan elimination.
    for i in 0..n {
        // Find the pivot; if it's too close to zero, the matrix is singular.
        let pivot = aug[[i, i]];
        if pivot.abs() < 1e-12 {
            return None;
        }
        // Normalize the pivot row.
        for j in 0..(2 * n) {
            aug[[i, j]] /= pivot;
        }
        // Eliminate the current column in all other rows.
        for k in 0..n {
            if k != i {
                let factor = aug[[k, i]];
                for j in 0..(2 * n) {
                    aug[[k, j]] -= factor * aug[[i, j]];
                }
            }
        }
    }

    // Extract the inverse from the augmented matrix.
    let mut inverse = Array2::<f64>::zeros((n, n));
    for i in 0..n {
        for j in 0..n {
            inverse[[i, j]] = aug[[i, n + j]];
        }
    }
    Some(inverse)
}

/// Radial Basis Function (RBF) kernel computes the similarity between two input vectors.
///
/// # Arguments
///
/// * `x1` - An Array1<f64> representing the first input vector.
/// * `x2` - An Array1<f64> representing the second input vector.
/// * `length_scale` - A f64 value that determines the smoothness of the kernel.
///
/// # Returns
///
/// A f64 value representing the kernel similarity.
fn rbf_kernel(x1: &Array1<f64>, x2: &Array1<f64>, length_scale: f64) -> f64 {
    let diff = x1 - x2;
    (-diff.dot(&diff) / (2.0 * length_scale.powi(2))).exp()
}

/// Constructs the kernel matrix for Gaussian Process regression using the RBF kernel.
///
/// # Arguments
///
/// * `train_data` - An Array2<f64> where each row is a training sample.
/// * `length_scale` - A f64 value for the kernel length scale.
///
/// # Returns
///
/// An Array2<f64> representing the kernel matrix.
fn build_kernel_matrix(train_data: &Array2<f64>, length_scale: f64) -> Array2<f64> {
    let n = train_data.nrows();
    let mut kernel_matrix = Array2::<f64>::zeros((n, n));
    for i in 0..n {
        for j in 0..n {
            // Convert the row views to owned arrays.
            let x1 = train_data.row(i).to_owned();
            let x2 = train_data.row(j).to_owned();
            kernel_matrix[[i, j]] = rbf_kernel(&x1, &x2, length_scale);
        }
    }
    kernel_matrix
}

/// Performs Gaussian Process prediction for new test data based on training data and outputs.
///
/// This function builds the kernel matrix for the training data, computes its inverse using our custom
/// inversion function, and then predicts the output for each test sample using the kernel vector between
/// the test sample and training samples.
///
/// # Arguments
///
/// * `train_data` - An Array2<f64> with training input samples.
/// * `train_outputs` - An Array1<f64> containing the training outputs.
/// * `test_data` - An Array2<f64> with test input samples.
/// * `length_scale` - A f64 value representing the kernel length scale.
///
/// # Returns
///
/// An Array1<f64> of predictions for the test data.
fn gp_predict(
    train_data: &Array2<f64>,
    train_outputs: &Array1<f64>,
    test_data: &Array2<f64>,
    length_scale: f64,
) -> Array1<f64> {
    // Build the kernel matrix for the training data.
    let kernel_matrix = build_kernel_matrix(train_data, length_scale);
    // Compute the inverse of the kernel matrix using our custom inversion routine.
    let inverse_kernel = invert_matrix(&kernel_matrix)
        .expect("Kernel matrix inversion failed; check if the matrix is well-conditioned");

    let mut predictions = Array1::<f64>::zeros(test_data.nrows());
    // Predict each test sample by computing the dot product between the kernel vector and the weights.
    for (i, test_point) in test_data.outer_iter().enumerate() {
        let mut k_star = Array1::<f64>::zeros(train_data.nrows());
        for (j, train_point) in train_data.outer_iter().enumerate() {
            k_star[j] = rbf_kernel(&train_point.to_owned(), &test_point.to_owned(), length_scale);
        }
        let prediction = k_star.dot(&inverse_kernel.dot(train_outputs));
        predictions[i] = prediction;
    }
    predictions
}

fn main() {
    // Simulate training data: for example, stress as a function of load.
    let num_train_points = 10;
    let load_values = Array1::linspace(0.0, 100.0, num_train_points);
    // Generate stress values using a linear relationship with added noise.
    // Note: The usage of `rand::thread_rng()` and `gen_range` here may trigger deprecation warnings.
    let mut rng = thread_rng();
    let stress_values = load_values.mapv(|load| 2.0 * load + 10.0 + rng.gen_range(-2.0..2.0));

    // Reshape the load values into a two-dimensional training data matrix.
    let train_data = load_values
        .clone()
        .into_shape((num_train_points, 1))
        .expect("Reshaping failed");
    let train_outputs = stress_values.clone();

    // Define test data: new load values where predictions are required.
    let num_test_points = 20;
    let test_load_values = Array1::linspace(0.0, 100.0, num_test_points);
    let test_data = test_load_values
        .clone()
        .into_shape((num_test_points, 1))
        .expect("Reshaping failed");

    // Specify the length scale for the RBF kernel.
    let length_scale = 10.0;
    // Use the Gaussian Process surrogate model to predict stress for the test data.
    let predicted_stress = gp_predict(&train_data, &train_outputs, &test_data, length_scale);

    // Display the predictions for each test load value.
    for (load, prediction) in test_load_values.iter().zip(predicted_stress.iter()) {
        println!("Load: {:.2} N, Predicted Stress: {:.2} Pa", load, prediction);
    }
}

use ndarray::Array1;
use rand::Rng;
use rand::thread_rng; // Note: may be deprecated in your version.
use rand_distr::{Normal, Distribution};

/// Performs a forward pass through a simple single-layer neural network that acts as a surrogate model.
/// 
/// In this model, the network consists of one weight vector and a bias term. The output is computed as the dot product
/// of the weights with the input plus the bias, followed by a Rectified Linear Unit (ReLU) activation to introduce nonlinearity.
/// 
/// # Arguments
/// 
/// * `weights` - An Array1<f64> representing the weight vector.
/// * `biases` - An Array1<f64> containing the bias; for simplicity, we use a single bias term.
/// * `inputs` - An Array1<f64> representing the input features (e.g., load).
/// 
/// # Returns
/// 
/// A f64 value representing the network's output (e.g., predicted stress).
fn neural_network_forward(weights: &Array1<f64>, biases: &Array1<f64>, inputs: &Array1<f64>) -> f64 {
    let weighted_sum = weights.dot(inputs) + biases[0];
    // Apply ReLU activation: returns the input if positive, else zero.
    weighted_sum.max(0.0)
}

fn main() {
    // Simulate training data for material behavior: stress as a function of load.
    let num_train_points = 10;
    let load_values = Array1::linspace(0.0, 100.0, num_train_points);
    
    // Generate corresponding stress values with a linear relationship plus some noise.
    let mut rng = thread_rng();
    let stress_values = load_values.mapv(|load| 2.0 * load + 10.0 + rng.gen_range(-2.0..2.0));
    // (stress_values is generated for demonstration; it's not used further in this example)

    // Prepare a normal distribution for generating weights and biases.
    let normal = Normal::new(0.0, 1.0).unwrap();
    
    // Manually generate a 1-element Array1<f64> for weights and biases.
    let weights: Array1<f64> = (0..1).map(|_| normal.sample(&mut rng)).collect();
    let biases: Array1<f64> = (0..1).map(|_| normal.sample(&mut rng)).collect();

    // Use the neural network surrogate model to predict stress for each training sample.
    for load in load_values.iter() {
        let input = Array1::from_elem(1, *load);
        let predicted_stress = neural_network_forward(&weights, &biases, &input);
        println!("Load: {:.2} N, Predicted Stress: {:.2} Pa", load, predicted_stress);
    }
}

use ndarray::Array1;
use rand::Rng;
use rand_distr::{Normal, Distribution};

/// Computes the structural response (stress) for a given fluid pressure and structural stiffness.
/// 
/// In this simplified model, the stress on the structure is assumed to be proportional to the fluid pressure
/// divided by the structural stiffness. This model serves as an approximation for how a structure might react to fluid forces.
/// 
/// # Arguments
///
/// * `pressure` - A f64 value representing the fluid pressure acting on the structure (in Pascals).
/// * `stiffness` - A f64 value representing the structural stiffness (in N/m).
///
/// # Returns
///
/// A f64 value representing the computed stress (in Pascals).
fn structural_model(pressure: f64, stiffness: f64) -> f64 {
    pressure / stiffness
}

/// Computes the fluid pressure using a simplified fluid dynamics model.
/// 
/// This function uses a basic form of the Bernoulli equation to compute the pressure based on fluid velocity and density.
/// The pressure is proportional to one-half times the density multiplied by the square of the velocity.
/// 
/// # Arguments
///
/// * `velocity` - A f64 value representing the fluid velocity (in m/s).
/// * `density` - A f64 value representing the fluid density (in kg/m³).
///
/// # Returns
///
/// A f64 value representing the fluid pressure (in Pascals).
fn fluid_model(velocity: f64, density: f64) -> f64 {
    0.5 * density * velocity.powi(2)
}

/// Performs Monte Carlo simulation to propagate uncertainties through a coupled fluid-structure model.
/// 
/// In this model, uncertainties in the fluid domain are introduced by assuming that fluid velocity and density are random variables
/// described by normal distributions with specified means and standard deviations. For each Monte Carlo sample, a random value is drawn
/// for velocity and density, the fluid pressure is computed using the fluid model, and then the structural response is computed using the
/// structural model. The result is a vector of structural responses that reflects the uncertainty propagated from the fluid properties.
/// 
/// # Arguments
///
/// * `velocity_mean` - The mean fluid velocity (in m/s).
/// * `velocity_std` - The standard deviation of the fluid velocity (in m/s).
/// * `density_mean` - The mean fluid density (in kg/m³).
/// * `density_std` - The standard deviation of the fluid density (in kg/m³).
/// * `stiffness` - The structural stiffness (in N/m).
/// * `num_samples` - The number of Monte Carlo samples to generate.
///
/// # Returns
///
/// A vector of f64 values representing the simulated structural responses (stresses).
fn monte_carlo_fsi(
    velocity_mean: f64,
    velocity_std: f64,
    density_mean: f64,
    density_std: f64,
    stiffness: f64,
    num_samples: usize
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    let velocity_dist = Normal::new(velocity_mean, velocity_std).unwrap();
    let density_dist = Normal::new(density_mean, density_std).unwrap();

    (0..num_samples)
        .map(|_| {
            // Sample a random fluid velocity and density from their respective distributions.
            let velocity_sample = velocity_dist.sample(&mut rng);
            let density_sample = density_dist.sample(&mut rng);

            // Compute the fluid pressure using the simplified fluid model.
            let pressure = fluid_model(velocity_sample, density_sample);

            // Compute the structural response (stress) based on the computed fluid pressure.
            structural_model(pressure, stiffness)
        })
        .collect()
}

/// Computes statistical measures (mean and standard deviation) of a set of simulation outputs.
/// 
/// This function calculates the mean value of the outputs and the standard deviation, providing a quantitative
/// measure of the uncertainty in the simulation results.
/// 
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the simulation outputs (e.g., structural responses).
///
/// # Returns
///
/// A tuple (mean, standard deviation) that summarizes the statistical distribution of the outputs.
fn compute_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / (samples.len() as f64);
    (mean, variance.sqrt())
}

fn main() {
    // Define the parameters for the fluid domain.
    let velocity_mean = 10.0;  // Mean fluid velocity in m/s.
    let velocity_std = 2.0;    // Standard deviation of fluid velocity.
    let density_mean = 1.0;    // Mean fluid density in kg/m³.
    let density_std = 0.1;     // Standard deviation of fluid density.

    // Define the parameter for the structural domain.
    let stiffness = 500.0;     // Structural stiffness in N/m.

    // Specify the number of Monte Carlo samples for uncertainty quantification.
    let num_samples = 1000;

    // Perform the Monte Carlo simulation for the fluid-structure interaction model.
    let structural_responses = monte_carlo_fsi(
        velocity_mean,
        velocity_std,
        density_mean,
        density_std,
        stiffness,
        num_samples
    );

    // Compute the mean and standard deviation of the structural responses to assess uncertainty.
    let (mean_response, std_response) = compute_statistics(&structural_responses);

    println!("Mean structural response: {:.2} N", mean_response);
    println!("Standard deviation of structural response: {:.2} N", std_response);
}

use rand::Rng;
use rand_distr::{Normal, Distribution};
use std::f64::consts::PI;

/// Defines a simple quantum model where the energy is modeled as a quadratic function of an unknown parameter.
///
/// In this example, the energy is computed as E = param² + 2 * param + 1, which serves as a simplified representation
/// of a quantum system's behavior. The model is used to predict energy given a value for the parameter.
///
/// # Arguments
///
/// * `param` - A f64 value representing the unknown parameter in the model.
///
/// # Returns
///
/// A f64 value representing the computed energy.
fn quantum_model(param: f64) -> f64 {
    param.powi(2) + 2.0 * param + 1.0
}

/// Computes the probability density function (pdf) of a normal (Gaussian) distribution.
///
/// # Arguments
///
/// * `x` - The value at which the pdf is evaluated.
/// * `mean` - The mean (μ) of the distribution.
/// * `sigma` - The standard deviation (σ) of the distribution.
///
/// # Returns
///
/// A f64 value representing the probability density at `x`.
fn normal_pdf(x: f64, mean: f64, sigma: f64) -> f64 {
    (1.0 / (sigma * (2.0 * PI).sqrt())) * (-0.5 * ((x - mean) / sigma).powi(2)).exp()
}

/// Computes the likelihood of observing the given data for a specific parameter value in the quantum model.
///
/// The likelihood is based on a Gaussian noise model with a fixed standard deviation, where the probability of
/// observing the data is evaluated using the normal pdf centered at the model prediction.
///
/// # Arguments
///
/// * `data` - A f64 value representing the observed energy measurement.
/// * `param` - A f64 value representing the parameter to be evaluated.
/// * `noise` - A f64 value representing the standard deviation of the measurement noise.
///
/// # Returns
///
/// A f64 value representing the likelihood of the data given the parameter.
fn likelihood(data: f64, param: f64, noise: f64) -> f64 {
    let prediction = quantum_model(param);
    normal_pdf(data, prediction, noise)
}

/// Uses the Metropolis-Hastings algorithm to sample from the posterior distribution of the parameter in the quantum model.
///
/// Starting from an initial guess (prior_mean), the algorithm proposes new parameter values and accepts them based on the
/// ratio of their likelihoods. This process is repeated for a specified number of iterations to generate a chain of samples
/// that approximate the posterior distribution.
///
/// # Arguments
///
/// * `data` - A f64 value representing the observed energy measurement.
/// * `prior_mean` - A f64 value representing the mean of the prior distribution for the parameter.
/// * `prior_std` - A f64 value representing the standard deviation of the prior distribution (used here to define the proposal range).
/// * `noise` - A f64 value representing the measurement noise standard deviation.
/// * `num_samples` - The number of MCMC samples to generate.
///
/// # Returns
///
/// A vector of f64 values representing samples from the posterior distribution of the parameter.
fn metropolis_hastings(data: f64, prior_mean: f64, prior_std: f64, noise: f64, num_samples: usize) -> Vec<f64> {
    // Note: `rand::thread_rng()` is marked deprecated; consider updating to the latest API if available.
    let mut rng = rand::thread_rng();
    let mut samples = Vec::with_capacity(num_samples);
    
    // Initialize the parameter with the prior mean.
    let mut current_param = prior_mean;

    // Use a Gaussian proposal distribution to generate new candidate parameters.
    let proposal_dist = Normal::new(0.0, prior_std).unwrap();

    // Iteratively sample the parameter space.
    for _ in 0..num_samples {
        // Propose a new parameter by adding a small perturbation.
        let proposed_param = current_param + proposal_dist.sample(&mut rng);

        // Calculate the likelihoods for the current and proposed parameter values.
        let current_likelihood = likelihood(data, current_param, noise);
        let proposed_likelihood = likelihood(data, proposed_param, noise);

        // Calculate the acceptance ratio.
        let acceptance_ratio = proposed_likelihood / current_likelihood;

        // Accept the proposed parameter with probability equal to the acceptance ratio.
        if rng.gen::<f64>() < acceptance_ratio {
            current_param = proposed_param;
        }

        // Record the current parameter value.
        samples.push(current_param);
    }

    samples
}

/// Computes statistical summaries (mean and standard deviation) of the posterior samples.
///
/// These statistics provide insight into the central tendency and variability of the estimated parameter,
/// allowing for an assessment of the uncertainty in the estimation process.
///
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the posterior samples.
///
/// # Returns
///
/// A tuple (mean, standard deviation) summarizing the posterior distribution.
fn compute_posterior_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / samples.len() as f64;
    (mean, variance.sqrt())
}

fn main() {
    // Simulated observation data: e.g., measured energy from a quantum system.
    let observed_data = 10.0;
    // Measurement noise standard deviation.
    let noise = 1.0;

    // Prior information for the parameter.
    let prior_mean = 2.0;
    let prior_std = 1.0;

    // Number of samples for the MCMC algorithm.
    let num_samples = 10_000;
    
    // Perform Metropolis-Hastings sampling to generate posterior samples.
    let samples = metropolis_hastings(observed_data, prior_mean, prior_std, noise, num_samples);

    // Compute the posterior mean and standard deviation from the samples.
    let (posterior_mean, posterior_std) = compute_posterior_statistics(&samples);

    // Output the results, including the estimated parameter and its uncertainty.
    println!("Posterior mean of the parameter: {:.4}", posterior_mean);
    println!("Posterior uncertainty (standard deviation): {:.4}", posterior_std);
}

use ndarray::Array1;

/// A simple structural model based on Hooke's law, which relates the displacement of a structure
/// to the applied force and the structural stiffness. This model is used to compute the displacement
/// given a force and a stiffness value.
///
/// # Arguments
///
/// * `force` - A f64 value representing the applied force in Newtons.
/// * `stiffness` - A f64 value representing the stiffness of the structure in N/m.
///
/// # Returns
///
/// A f64 value representing the computed displacement in meters.
fn structural_model(force: f64, stiffness: f64) -> f64 {
    force / stiffness
}

/// Performs solution verification by examining the convergence behavior of the computed displacement
/// as the force step size is refined. This function calculates displacements for a set of force values,
/// then refines the force increments and computes displacements again to assess if the solution converges.
///
/// # Arguments
///
/// * `stiffness` - A f64 value representing the structural stiffness (in N/m).
fn verify_convergence(stiffness: f64) {
    // Generate a vector of force values using a coarse step size.
    let force_values: Vec<f64> = (1..=10).map(|i| i as f64 * 10.0).collect();
    let mut displacements = Array1::<f64>::zeros(force_values.len());

    // Compute displacements for each force value.
    for (i, &force) in force_values.iter().enumerate() {
        displacements[i] = structural_model(force, stiffness);
    }

    // Output the displacements for the coarse grid.
    println!("Displacements for coarse force steps:");
    for (force, displacement) in force_values.iter().zip(displacements.iter()) {
        println!("Force: {:.2} N, Displacement: {:.4} m", force, displacement);
    }

    // Refine the force step size by halving it.
    let refined_force_values: Vec<f64> = (1..=20).map(|i| i as f64 * 5.0).collect();
    let mut refined_displacements = Array1::<f64>::zeros(refined_force_values.len());

    // Compute displacements for the refined force values.
    for (i, &force) in refined_force_values.iter().enumerate() {
        refined_displacements[i] = structural_model(force, stiffness);
    }

    // Output the refined displacements to check for convergence.
    println!("\nDisplacements for refined force steps:");
    for (force, displacement) in refined_force_values.iter().zip(refined_displacements.iter()) {
        println!("Force: {:.2} N, Displacement: {:.4} m", force, displacement);
    }
}

fn main() {
    let stiffness = 100.0; // Structural stiffness in N/m

    // Execute solution verification by analyzing convergence with different force step sizes.
    verify_convergence(stiffness);
}

use rand::Rng;

/// A simple thermal model that predicts the temperature at a given distance from a heat source.
/// In this model, temperature decreases with the square of the distance, modulated by a thermal conductivity parameter.
///
/// # Arguments
///
/// * `distance` - A f64 value representing the distance from the heat source (in meters).
/// * `thermal_conductivity` - A f64 value representing the thermal conductivity (in appropriate units).
///
/// # Returns
///
/// A f64 value representing the predicted temperature in Celsius.
fn thermal_model(distance: f64, thermal_conductivity: f64) -> f64 {
    thermal_conductivity / (1.0 + distance.powi(2))
}

/// Validates the thermal model by comparing its predictions with experimental data.
/// This function takes a list of experimental measurements (distance, temperature pairs) and computes the
/// predicted temperature for each distance using the thermal model. It then prints a side-by-side comparison
/// of simulated and experimental temperatures.
///
/// # Arguments
///
/// * `experimental_data` - A slice of tuples where each tuple contains (distance, experimental temperature).
/// * `thermal_conductivity` - A f64 value representing the thermal conductivity used in the simulation.
fn validate_against_experiment(experimental_data: &[(f64, f64)], thermal_conductivity: f64) {
    println!("Comparing simulation results to experimental data:");
    for &(distance, experimental_temp) in experimental_data.iter() {
        let simulated_temp = thermal_model(distance, thermal_conductivity);
        println!(
            "Distance: {:.2} m, Simulated Temp: {:.2} °C, Experimental Temp: {:.2} °C",
            distance, simulated_temp, experimental_temp
        );
    }
}

fn main() {
    // Define a set of experimental data points (distance in meters, temperature in Celsius).
    let experimental_data = vec![
        (0.0, 100.0),
        (0.5, 60.0),
        (1.0, 30.0),
        (1.5, 20.0),
        (2.0, 15.0),
    ];

    // Define the thermal conductivity parameter used in the simulation.
    let thermal_conductivity = 80.0;

    // Validate the thermal model by comparing simulated results with the experimental data.
    validate_against_experiment(&experimental_data, thermal_conductivity);
}

use rand::Rng;
use rand_distr::{Normal, Distribution};
use ndarray::Array1;

/// A simplified model for reactor performance that relates the reactor's performance to thermal conductivity and degradation.
/// In this model, reactor performance decreases with lower thermal conductivity and increases with lower degradation.
/// The model is formulated as performance = thermal_conductivity / (1.0 + degradation).
///
/// # Arguments
///
/// * `thermal_conductivity` - A f64 value representing the thermal conductivity (in W/mK).
/// * `degradation` - A f64 value representing the degradation factor (unitless).
///
/// # Returns
///
/// A f64 value representing the reactor performance.
fn reactor_performance(thermal_conductivity: f64, degradation: f64) -> f64 {
    thermal_conductivity / (1.0 + degradation)
}

/// Performs Monte Carlo simulation to propagate uncertainties in thermal conductivity and degradation through the reactor performance model.
///
/// This function samples thermal conductivity and degradation from their respective normal distributions,
/// evaluates the reactor performance for each sample, and returns a vector of performance outcomes.
///
/// # Arguments
///
/// * `thermal_mean` - Mean thermal conductivity (in W/mK).
/// * `thermal_std` - Standard deviation of thermal conductivity.
/// * `degradation_mean` - Mean degradation factor.
/// * `degradation_std` - Standard deviation of degradation.
/// * `num_samples` - Number of Monte Carlo samples to generate.
///
/// # Returns
///
/// A vector of f64 values representing the simulated reactor performance.
fn monte_carlo_uq(
    thermal_mean: f64,
    thermal_std: f64,
    degradation_mean: f64,
    degradation_std: f64,
    num_samples: usize,
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    let thermal_dist = Normal::new(thermal_mean, thermal_std).unwrap();
    let degradation_dist = Normal::new(degradation_mean, degradation_std).unwrap();

    (0..num_samples)
        .map(|_| {
            // Sample uncertain parameters from their distributions.
            let thermal_sample = thermal_dist.sample(&mut rng);
            let degradation_sample = degradation_dist.sample(&mut rng);
            // Compute reactor performance using the sampled parameters.
            reactor_performance(thermal_sample, degradation_sample)
        })
        .collect()
}

/// Computes statistical measures (mean and standard deviation) of the performance samples.
///
/// This function calculates the average performance and its standard deviation, providing a measure
/// of the central tendency and spread of the outcomes, which reflects the propagated uncertainty.
///
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the performance outcomes from Monte Carlo simulation.
///
/// # Returns
///
/// A tuple (mean, standard deviation) summarizing the performance distribution.
fn compute_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / samples.len() as f64;
    (mean, variance.sqrt())
}

fn main() {
    // Define uncertain parameters for the reactor model.
    let thermal_mean = 50.0;   // Mean thermal conductivity (W/mK).
    let thermal_std = 5.0;     // Standard deviation of thermal conductivity.
    let degradation_mean = 0.1; // Mean degradation factor.
    let degradation_std = 0.05; // Standard deviation of degradation.

    // Specify the number of Monte Carlo samples.
    let num_samples = 10000;

    // Run the Monte Carlo simulation to obtain performance samples.
    let performance_samples = monte_carlo_uq(
        thermal_mean,
        thermal_std,
        degradation_mean,
        degradation_std,
        num_samples,
    );

    // Compute the mean and uncertainty (standard deviation) of the reactor performance.
    let (mean_performance, std_performance) = compute_statistics(&performance_samples);

    println!("Mean reactor performance: {:.4}", mean_performance);
    println!("Uncertainty (standard deviation): {:.4}", std_performance);
}

use rand::Rng;
use rand_distr::{Normal, Distribution};

/// A simplified climate model that estimates the temperature increase based on contributions from atmospheric and oceanic parameters.
/// The model weights the atmospheric contribution at 70% and the oceanic contribution at 30% to compute the overall temperature increase.
///
/// # Arguments
///
/// * `atmosphere` - A f64 value representing an atmospheric parameter.
/// * `ocean` - A f64 value representing an oceanic parameter.
///
/// # Returns
///
/// A f64 value representing the predicted temperature increase (in °C).
fn climate_model(atmosphere: f64, ocean: f64) -> f64 {
    0.7 * atmosphere + 0.3 * ocean
}

/// Performs ensemble simulations to propagate uncertainties in atmospheric and oceanic parameters through the climate model.
///
/// This function generates multiple ensembles by sampling atmospheric and oceanic parameters from their respective normal distributions,
/// then computes the predicted temperature increase for each ensemble. The ensemble of results reflects the uncertainty in the predictions.
///
/// # Arguments
///
/// * `atmos_mean` - Mean atmospheric parameter.
/// * `atmos_std` - Standard deviation of the atmospheric parameter.
/// * `ocean_mean` - Mean oceanic parameter.
/// * `ocean_std` - Standard deviation of the oceanic parameter.
/// * `num_ensembles` - Number of ensemble simulations to perform.
///
/// # Returns
///
/// A vector of f64 values representing the ensemble predictions of temperature increase.
fn ensemble_uq(
    atmos_mean: f64,
    atmos_std: f64,
    ocean_mean: f64,
    ocean_std: f64,
    num_ensembles: usize,
) -> Vec<f64> {
    let mut rng = rand::thread_rng();
    let atmos_dist = Normal::new(atmos_mean, atmos_std).unwrap();
    let ocean_dist = Normal::new(ocean_mean, ocean_std).unwrap();

    (0..num_ensembles)
        .map(|_| {
            // Sample atmospheric and oceanic parameters.
            let atmosphere_sample = atmos_dist.sample(&mut rng);
            let ocean_sample = ocean_dist.sample(&mut rng);
            // Compute the predicted temperature increase using the climate model.
            climate_model(atmosphere_sample, ocean_sample)
        })
        .collect()
}

/// Computes the mean and standard deviation of ensemble simulation outputs.
///
/// This function aggregates the ensemble predictions to determine the average predicted temperature increase
/// and the uncertainty associated with the prediction.
///
/// # Arguments
///
/// * `samples` - A slice of f64 values representing the ensemble outputs.
///
/// # Returns
///
/// A tuple (mean, standard deviation) summarizing the ensemble prediction.
fn compute_statistics(samples: &[f64]) -> (f64, f64) {
    let mean: f64 = samples.iter().sum::<f64>() / samples.len() as f64;
    let variance: f64 = samples.iter().map(|&x| (x - mean).powi(2)).sum::<f64>() / samples.len() as f64;
    (mean, variance.sqrt())
}

fn main() {
    // Define the uncertain atmospheric and oceanic parameters.
    let atmos_mean = 1.0;   // Mean atmospheric parameter.
    let atmos_std = 0.2;    // Standard deviation of the atmospheric parameter.
    let ocean_mean = 0.8;   // Mean oceanic parameter.
    let ocean_std = 0.1;    // Standard deviation of the oceanic parameter.

    // Specify the number of ensemble simulations.
    let num_ensembles = 5000;

    // Run ensemble simulations to propagate uncertainties through the climate model.
    let temperature_samples = ensemble_uq(atmos_mean, atmos_std, ocean_mean, ocean_std, num_ensembles);

    // Compute the mean predicted temperature increase and its uncertainty.
    let (mean_temperature, std_temperature) = compute_statistics(&temperature_samples);

    println!("Mean predicted temperature increase: {:.4} °C", mean_temperature);
    println!("Uncertainty (standard deviation): {:.4} °C", std_temperature);
}