Struct Uncertain 

pub struct Uncertain<S>
where
    S: Scalar,
{
    pub mean: S,
    pub variance: S,
}

An uncertain value with mean and variance.

Example

use numra_core::uncertainty::Uncertain;

// A measurement of 10.0 ± 0.5 (std dev)
let x: Uncertain<f64> = Uncertain::new(10.0, 0.25); // variance = 0.5²

// Scale by 2: mean scales, variance scales by 4
let y = x.scale(2.0);
assert!((y.mean - 20.0).abs() < 1e-10);
assert!((y.variance - 1.0).abs() < 1e-10);

Fields

mean: S

Mean value

variance: S

Variance (σ²)

Implementations

impl<S> Uncertain<S>

where S: Scalar,

pub fn new(mean: S, variance: S) -> Uncertain<S>

Create a new uncertain value.

pub fn from_std(mean: S, std: S) -> Uncertain<S>

Create from mean and standard deviation.

Examples found in repository

examples/lorenz.rs (line 46)

fn main() {
    println!("=== Lorenz Attractor with Uncertainty ===\n");

    // Standard Lorenz parameters
    let params = LorenzParams::default();

    // Initial conditions (with uncertainty)
    let x0 = Uncertain::from_std(1.0, 0.1); // x = 1.0 ± 0.1
    let y0 = Uncertain::from_std(1.0, 0.1); // y = 1.0 ± 0.1
    let z0 = Uncertain::from_std(1.0, 0.1); // z = 1.0 ± 0.1

    println!("Initial conditions:");
    println!("  x₀ = {:.3} ± {:.3}", x0.mean, x0.std());
    println!("  y₀ = {:.3} ± {:.3}", y0.mean, y0.std());
    println!("  z₀ = {:.3} ± {:.3}", z0.mean, z0.std());
    println!();

    // Create the ODE problem
    let problem = OdeProblem::new(
        |_t, y: &[f64], dydt: &mut [f64]| {
            let (x, y_val, z) = (y[0], y[1], y[2]);
            let sigma = 10.0;
            let rho = 28.0;
            let beta = 8.0 / 3.0;

            dydt[0] = sigma * (y_val - x);
            dydt[1] = x * (rho - z) - y_val;
            dydt[2] = x * y_val - beta * z;
        },
        0.0,
        20.0,
        vec![x0.mean, y0.mean, z0.mean],
    );

    // Solve with default options
    let options = SolverOptions::default().rtol(1e-8).atol(1e-10);

    println!("Solving Lorenz system from t=0 to t=20...");
    let result = DoPri5::solve(&problem, 0.0, 20.0, &[x0.mean, y0.mean, z0.mean], &options)
        .expect("Solver failed");

    // Get solution statistics
    println!("  Steps accepted: {}", result.stats.n_accept);
    println!("  Steps rejected: {}", result.stats.n_reject);
    println!("  Function evals: {}", result.stats.n_eval);
    println!();

    // Final state
    let y_final = result.y_final().expect("No final state");
    println!("Final state at t=20:");
    println!("  x = {:.6}", y_final[0]);
    println!("  y = {:.6}", y_final[1]);
    println!("  z = {:.6}", y_final[2]);
    println!();

    // Trajectory summary
    println!("Trajectory has {} points", result.t.len());
    println!();

    // === Sensitivity Analysis ===
    println!("=== Sensitivity Analysis ===\n");

    // Compute sensitivity of final z-coordinate with respect to parameters
    let solve_for_z = |p: &[f64]| -> f64 {
        let problem = OdeProblem::new(
            move |_t, y: &[f64], dydt: &mut [f64]| {
                let (x, y_val, z) = (y[0], y[1], y[2]);
                dydt[0] = p[0] * (y_val - x); // σ
                dydt[1] = x * (p[1] - z) - y_val; // ρ
                dydt[2] = x * y_val - p[2] * z; // β
            },
            0.0,
            5.0, // Shorter time for sensitivity (chaos makes long-time sensitivity ill-defined)
            vec![1.0, 1.0, 1.0],
        );

        let options = SolverOptions::default().rtol(1e-6);
        let result = DoPri5::solve(&problem, 0.0, 5.0, &[1.0, 1.0, 1.0], &options)
            .expect("Solver failed in sensitivity");
        result.y_final().unwrap()[2] // Return final z
    };

    let nominal_params = [params.sigma, params.rho, params.beta];
    let param_names = ["σ", "ρ", "β"];

    let sens_result = compute_sensitivities(solve_for_z, &nominal_params, &param_names, None);

    println!("Sensitivity of z(t=5) to parameters:");
    for s in &sens_result.sensitivities {
        println!(
            "  ∂z/∂{}: {:.4} (normalized: {:.4})",
            s.name, s.coefficient, s.normalized
        );
    }
    println!();

    // Find most sensitive parameter
    if let Some(most_sens) = sens_result.most_sensitive() {
        println!(
            "Most sensitive parameter: {} (|normalized| = {:.4})",
            most_sens.name,
            most_sens.normalized.abs()
        );
    }
    println!();

    // === Uncertainty Propagation ===
    println!("=== Uncertainty Propagation ===\n");

    // Assume 5% uncertainty on each parameter
    let param_std = [0.5, 1.4, 0.133]; // 5% of nominal
    let param_vars: Vec<f64> = param_std.iter().map(|s| s * s).collect();

    let total_var = sens_result.propagate_uncertainty(&param_vars);
    let total_std = total_var.sqrt();

    println!("Parameter uncertainties:");
    println!("  σ = {:.2} ± {:.3}", params.sigma, param_std[0]);
    println!("  ρ = {:.2} ± {:.3}", params.rho, param_std[1]);
    println!("  β = {:.4} ± {:.4}", params.beta, param_std[2]);
    println!();

    println!("Propagated uncertainty in z(t=5):");
    println!("  z = {:.4} ± {:.4}", sens_result.output, total_std);
    println!();

    // Individual contributions
    println!("Uncertainty contributions:");
    for (i, s) in sens_result.sensitivities.iter().enumerate() {
        let contribution = (s.coefficient * param_std[i]).abs();
        let pct = 100.0 * contribution * contribution / total_var;
        println!("  {}: {:.2}%", s.name, pct);
    }
    println!();

    // === Ensemble of trajectories ===
    println!("=== Ensemble with Perturbed Initial Conditions ===\n");

    let n_samples = 5;
    println!("Running {} trajectories with perturbed ICs...", n_samples);

    let perturbations = [
        [1.0, 1.0, 1.0],
        [1.1, 1.0, 1.0],
        [1.0, 1.1, 1.0],
        [0.9, 1.0, 1.0],
        [1.0, 0.9, 1.0],
    ];

    for ic in perturbations.iter() {
        let result = DoPri5::solve(&problem, 0.0, 10.0, ic, &options).expect("Solver failed");
        let y_end = result.y_final().unwrap();
        println!(
            "  IC [{:.1}, {:.1}, {:.1}] → z(10) = {:.4}",
            ic[0], ic[1], ic[2], y_end[2]
        );
    }
    println!();

    println!("Note: The Lorenz system is chaotic - small perturbations in ICs");
    println!("lead to exponentially diverging trajectories over time.");
    println!();

    println!("=== Done ===");
}

pub fn certain(value: S) -> Uncertain<S>

Create a certain value (zero variance).

pub fn std(&self) -> S

Standard deviation.

Examples found in repository

examples/lorenz.rs (line 51)

fn main() {
    println!("=== Lorenz Attractor with Uncertainty ===\n");

    // Standard Lorenz parameters
    let params = LorenzParams::default();

    // Initial conditions (with uncertainty)
    let x0 = Uncertain::from_std(1.0, 0.1); // x = 1.0 ± 0.1
    let y0 = Uncertain::from_std(1.0, 0.1); // y = 1.0 ± 0.1
    let z0 = Uncertain::from_std(1.0, 0.1); // z = 1.0 ± 0.1

    println!("Initial conditions:");
    println!("  x₀ = {:.3} ± {:.3}", x0.mean, x0.std());
    println!("  y₀ = {:.3} ± {:.3}", y0.mean, y0.std());
    println!("  z₀ = {:.3} ± {:.3}", z0.mean, z0.std());
    println!();

    // Create the ODE problem
    let problem = OdeProblem::new(
        |_t, y: &[f64], dydt: &mut [f64]| {
            let (x, y_val, z) = (y[0], y[1], y[2]);
            let sigma = 10.0;
            let rho = 28.0;
            let beta = 8.0 / 3.0;

            dydt[0] = sigma * (y_val - x);
            dydt[1] = x * (rho - z) - y_val;
            dydt[2] = x * y_val - beta * z;
        },
        0.0,
        20.0,
        vec![x0.mean, y0.mean, z0.mean],
    );

    // Solve with default options
    let options = SolverOptions::default().rtol(1e-8).atol(1e-10);

    println!("Solving Lorenz system from t=0 to t=20...");
    let result = DoPri5::solve(&problem, 0.0, 20.0, &[x0.mean, y0.mean, z0.mean], &options)
        .expect("Solver failed");

    // Get solution statistics
    println!("  Steps accepted: {}", result.stats.n_accept);
    println!("  Steps rejected: {}", result.stats.n_reject);
    println!("  Function evals: {}", result.stats.n_eval);
    println!();

    // Final state
    let y_final = result.y_final().expect("No final state");
    println!("Final state at t=20:");
    println!("  x = {:.6}", y_final[0]);
    println!("  y = {:.6}", y_final[1]);
    println!("  z = {:.6}", y_final[2]);
    println!();

    // Trajectory summary
    println!("Trajectory has {} points", result.t.len());
    println!();

    // === Sensitivity Analysis ===
    println!("=== Sensitivity Analysis ===\n");

    // Compute sensitivity of final z-coordinate with respect to parameters
    let solve_for_z = |p: &[f64]| -> f64 {
        let problem = OdeProblem::new(
            move |_t, y: &[f64], dydt: &mut [f64]| {
                let (x, y_val, z) = (y[0], y[1], y[2]);
                dydt[0] = p[0] * (y_val - x); // σ
                dydt[1] = x * (p[1] - z) - y_val; // ρ
                dydt[2] = x * y_val - p[2] * z; // β
            },
            0.0,
            5.0, // Shorter time for sensitivity (chaos makes long-time sensitivity ill-defined)
            vec![1.0, 1.0, 1.0],
        );

        let options = SolverOptions::default().rtol(1e-6);
        let result = DoPri5::solve(&problem, 0.0, 5.0, &[1.0, 1.0, 1.0], &options)
            .expect("Solver failed in sensitivity");
        result.y_final().unwrap()[2] // Return final z
    };

    let nominal_params = [params.sigma, params.rho, params.beta];
    let param_names = ["σ", "ρ", "β"];

    let sens_result = compute_sensitivities(solve_for_z, &nominal_params, &param_names, None);

    println!("Sensitivity of z(t=5) to parameters:");
    for s in &sens_result.sensitivities {
        println!(
            "  ∂z/∂{}: {:.4} (normalized: {:.4})",
            s.name, s.coefficient, s.normalized
        );
    }
    println!();

    // Find most sensitive parameter
    if let Some(most_sens) = sens_result.most_sensitive() {
        println!(
            "Most sensitive parameter: {} (|normalized| = {:.4})",
            most_sens.name,
            most_sens.normalized.abs()
        );
    }
    println!();

    // === Uncertainty Propagation ===
    println!("=== Uncertainty Propagation ===\n");

    // Assume 5% uncertainty on each parameter
    let param_std = [0.5, 1.4, 0.133]; // 5% of nominal
    let param_vars: Vec<f64> = param_std.iter().map(|s| s * s).collect();

    let total_var = sens_result.propagate_uncertainty(&param_vars);
    let total_std = total_var.sqrt();

    println!("Parameter uncertainties:");
    println!("  σ = {:.2} ± {:.3}", params.sigma, param_std[0]);
    println!("  ρ = {:.2} ± {:.3}", params.rho, param_std[1]);
    println!("  β = {:.4} ± {:.4}", params.beta, param_std[2]);
    println!();

    println!("Propagated uncertainty in z(t=5):");
    println!("  z = {:.4} ± {:.4}", sens_result.output, total_std);
    println!();

    // Individual contributions
    println!("Uncertainty contributions:");
    for (i, s) in sens_result.sensitivities.iter().enumerate() {
        let contribution = (s.coefficient * param_std[i]).abs();
        let pct = 100.0 * contribution * contribution / total_var;
        println!("  {}: {:.2}%", s.name, pct);
    }
    println!();

    // === Ensemble of trajectories ===
    println!("=== Ensemble with Perturbed Initial Conditions ===\n");

    let n_samples = 5;
    println!("Running {} trajectories with perturbed ICs...", n_samples);

    let perturbations = [
        [1.0, 1.0, 1.0],
        [1.1, 1.0, 1.0],
        [1.0, 1.1, 1.0],
        [0.9, 1.0, 1.0],
        [1.0, 0.9, 1.0],
    ];

    for ic in perturbations.iter() {
        let result = DoPri5::solve(&problem, 0.0, 10.0, ic, &options).expect("Solver failed");
        let y_end = result.y_final().unwrap();
        println!(
            "  IC [{:.1}, {:.1}, {:.1}] → z(10) = {:.4}",
            ic[0], ic[1], ic[2], y_end[2]
        );
    }
    println!();

    println!("Note: The Lorenz system is chaotic - small perturbations in ICs");
    println!("lead to exponentially diverging trajectories over time.");
    println!();

    println!("=== Done ===");
}

pub fn cv(&self) -> S

Coefficient of variation (std/mean).

pub fn scale(&self, a: S) -> Uncertain<S>

Scale by a constant: y = a*x. Variance scales by a².

pub fn add_const(&self, c: S) -> Uncertain<S>

Add a constant: y = x + c. Variance unchanged.

pub fn add(&self, other: &Uncertain<S>) -> Uncertain<S>

Add another uncertain value (assuming independence). y = x1 + x2, Var(y) = Var(x1) + Var(x2).

pub fn sub(&self, other: &Uncertain<S>) -> Uncertain<S>

Subtract another uncertain value (assuming independence).

pub fn mul(&self, other: &Uncertain<S>) -> Uncertain<S>

Multiply by another uncertain value (first-order approximation). For y = x1 * x2: E[y] ≈ E[x1] * E[x2] Var(y) ≈ E[x2]² * Var(x1) + E[x1]² * Var(x2)

pub fn div(&self, other: &Uncertain<S>) -> Uncertain<S>

Divide by another uncertain value (first-order approximation). For y = x1 / x2: E[y] ≈ E[x1] / E[x2] Var(y) ≈ (1/E[x2]²) * Var(x1) + (E[x1]/E[x2]²)² * Var(x2)

pub fn apply<F, D>(&self, f: F, df: D) -> Uncertain<S>

where F: Fn(S) -> S, D: Fn(S) -> S,

Apply a function with known derivative. For y = f(x): Var(y) ≈ (f’(x))² * Var(x).

pub fn square(&self) -> Uncertain<S>

Square: y = x².

pub fn sqrt_unc(&self) -> Uncertain<S>

Square root: y = √x.

pub fn exp(&self) -> Uncertain<S>

Exponential: y = e^x.

pub fn ln(&self) -> Uncertain<S>

Natural log: y = ln(x).

pub fn sin(&self) -> Uncertain<S>

Sine: y = sin(x).

pub fn cos(&self) -> Uncertain<S>

Cosine: y = cos(x).

Trait Implementations

impl<S> Clone for Uncertain<S>

where S: Clone + Scalar,

fn clone(&self) -> Uncertain<S>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<S> Debug for Uncertain<S>

where S: Debug + Scalar,

fn fmt(&self, f: &mut Formatter<'_>) -> Result<(), Error>

Formats the value using the given formatter.

impl<S> Default for Uncertain<S>

where S: Scalar,

fn default() -> Uncertain<S>

Returns the “default value” for a type.

impl<S> PartialEq for Uncertain<S>

where S: PartialEq + Scalar,

fn eq(&self, other: &Uncertain<S>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl<S> Copy for Uncertain<S>

where S: Copy + Scalar,

impl<S> StructuralPartialEq for Uncertain<S>

where S: Scalar,