// Copyright 2019 The Fuchsia Authors. All rights reserved.
// Use of this source code is governed by a BSD-style license that can be
// found in the LICENSE file.

use anyhow::{format_err, Error};

/// Runge-Kutta Fehlberg 4(5) parameters
///
/// These are readily available in literature, e.g.
/// <https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method>.
mod rkf45_params {
    /// Number of stages
    pub const NUM_STAGES: usize = 6;

    /// Runge-Kutta matrix
    pub static A: [[f64; NUM_STAGES - 1]; NUM_STAGES - 1] = [
        [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0],
        [3.0 / 32.0, 9.0 / 32.0, 0.0, 0.0, 0.0],
        [1932.0 / 2197.0, -7200.0 / 2197.0, 7296.0 / 2197.0, 0.0, 0.0],
        [439.0 / 216.0, -8.0, 3680.0 / 513.0, -845.0 / 4104.0, 0.0],
        [-8.0 / 27.0, 2.0, -3544.0 / 2565.0, 1859.0 / 4104.0, -11.0 / 40.0],
    ];

    /// 4th-order weights
    pub static B4: [f64; NUM_STAGES] =
        [25.0 / 216.0, 0.0, 1408.0 / 2565.0, 2197.0 / 4104.0, -1.0 / 5.0, 0.0];

    /// 5th-order weights
    pub static B5: [f64; NUM_STAGES] =
        [16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0];

    /// Nodes
    pub static C: [f64; NUM_STAGES - 1] = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0];
}

/// Parameters for adaptive time-stepping logic.
///
/// These incorporate a variety of heuristics and are mostly described in Press and Teukolsky's
/// [Adaptive Stepsize Runge-Kutta Integration, Computers in Physics 6, 188
/// (1992)](https://doi.org/10.1063/1.4823060).
mod adaptive_stepping_params {
    /// Safety factor to use when modifying time steps. This makes a revised time step slightly
    /// smaller than required by formal analysis, providing margin for inaccuracy in the error
    /// estimate.
    pub static SAFETY_FACTOR: f64 = 0.9;

    /// Only refine/coarsen the time step when the ratio of estimate error to desired error
    /// crosses these thresholds. Doing so prevents excessive recalculation of time steps.
    pub static ERROR_RATIO_REFINING_THRESHOLD: f64 = 1.1;
    pub static ERROR_RATIO_COARSENING_THRESHOLD: f64 = 0.5;

    /// RKF45's local truncation (single-step) error can be written as
    ///     E(dt) = C*dt^5 + O(dt^6).
    ///
    /// Suppose the desired error D is fixed with respect to dt. Then the optimal step size
    /// approximately satisfies
    ///     C*dt_opt^5 = D  ==> dt_opt = (D / C)^(1/5).
    /// If we have attempted a time step of size dt, then we can approximate
    ///     C = E(dt) / dt^5.
    /// This yields
    ///     dt_opt = dt * (D / E(dt))^(1/5) = (E(dt) / D)^(-1/5).
    ///
    /// Suppose instead that D is proprotional to dt, so D(dt) = D1*dt. Then the optimal step size
    /// satisfies
    ///     C*dt_opt^5 = D1*dt_opt  ==>  C = (D1 / C)^(1/4).
    /// If we have attempted a time step of size dt, with estimated error E(dt) and desired error
    /// D(dt), then we observe that
    ///     D(dt) / E(dt) = (D1*dt) / (C*dt^5)
    ///     ==> D1 / C = dt^4 (D(dt) / E(dt)).
    /// Using this relationship above yields
    ///     dt_opt = dt * (D(dt) / E(dt))^(1/4) = (E(dt) / D(dt))^(-1/4).
    ///
    /// Following the guidance of Press and Teukolsky, rather than inspect the form of the desired
    /// error, we simply use the more conservative exponent depending on the situation. When
    /// refining, we will have E/D > 1. Raising to the smaller exponent (-0.25) will yield a
    /// smaller time step. Converseley, when coarsening we have E/D < 1, and the larger exponent
    /// (-0.20) will yield a smaller timestep.
    pub static REFINING_EXPONENT: f64 = -0.25;
    pub static COARSENING_EXPONENT: f64 = -0.20;

    /// Bounds for the factor used when refining/coarsening the time step, to avoid too much of
    /// a change at once.
    pub static MIN_REFINING_FACTOR: f64 = 0.2;
    pub static MAX_COARSENING_FACTOR: f64 = 5.0;
}

/// Vector operations to implement on `[f64]` for convenience.
///
/// All methods involving another [f64] expect the two slices to be of equal length. The caller
/// is responsible for ensuring this. (Practically speaking, the slices involved are from `Vec`s
/// created using the length of `rkf45_adaptive`'s input `y`.)
trait VectorOperations {
    /// Scale by a scalar.
    fn scale(&mut self, a: f64);

    /// Copy to this vector from another one. Slice lengths must be equal.
    fn copy_from(&mut self, x: &Self);

    /// Add a vector to this one. Slice lengths must be equal.
    fn add(&mut self, x: &Self);

    /// Add a scalar times another vector. Slice lengths must be equal.
    fn add_ax(&mut self, a: f64, x: &Self);

    /// Subtract a vector from this one. Slice lengths must be equal.
    fn subtract(&mut self, x: &Self);

    /// Sets all elements to zero.
    fn to_zeros(&mut self);
}

impl VectorOperations for [f64] {
    fn scale(&mut self, a: f64) {
        self.iter_mut().for_each(|p| *p = *p * a);
    }

    fn copy_from(&mut self, x: &Self) {
        self.iter_mut().zip(x.iter()).for_each(|(p, q)| *p = *q);
    }

    fn add(&mut self, x: &Self) {
        self.iter_mut().zip(x.iter()).for_each(|(p, q)| *p += *q);
    }

    fn add_ax(&mut self, a: f64, x: &Self) {
        self.iter_mut().zip(x.iter()).for_each(|(p, q)| *p += a * *q);
    }

    fn subtract(&mut self, x: &Self) {
        self.iter_mut().zip(x.iter()).for_each(|(p, q)| *p -= *q);
    }

    fn to_zeros(&mut self) {
        self.iter_mut().for_each(|p| *p = 0.0)
    }
}

/// Takes a single time step using RKF45.
///
/// Args:
/// - `y`: Slice containing solution values at the beginning of the time step. This will be
///    updated in-place.
/// - `dydt`: Function that evaluates the time derivative of `y`.
/// - `tn`: Time value at the beginning of the step.
/// - `dt`: Length of the step.
///
/// Returns:
/// - Vector estimating elementwise numerical error over this single step, and the maximum ratio of
///   estimated error to error bound.
fn rkf45_step(
    y: &mut [f64],
    dydt: &impl Fn(f64, &[f64]) -> Vec<f64>,
    tn: f64,
    dt: f64,
    error_control: &ErrorControlOptions,
) -> (Vec<f64>, f64) {
    // Work array reused by several loops below.
    let mut work = vec![0.0; y.len()];

    let mut k = Vec::with_capacity(6);
    k.push(dydt(tn, y));

    use rkf45_params as params;
    for i in 1..params::NUM_STAGES {
        work.to_zeros();
        for j in 0..i {
            work.add_ax(params::A[i - 1][j], &k[j]);
        }
        // Effectively, k[i] = dydt(tn + C[i-1] * dt, y + dt * work).
        work.scale(dt);
        work.add(&y);
        k.push(dydt(tn + params::C[i - 1] * dt, &work));
    }

    // The 5th-order solution will be used for error-estimation.
    let mut y_5th_order = y.to_vec();
    work.to_zeros();
    for i in 0..params::B5.len() {
        work.add_ax(params::B5[i], &k[i]);
    }
    y_5th_order.add_ax(dt, &work);

    // y is updated with the 4th-order solution.
    work.to_zeros();
    for i in 0..params::B4.len() {
        work.add_ax(params::B4[i], &k[i]);
    }
    y.add_ax(dt, &work);

    let mut error_estimate = y_5th_order;
    error_estimate.subtract(&y);
    error_estimate.iter_mut().for_each(|x| *x = (*x).abs());

    let mut max_error_ratio = 0.0;
    for i in 0..y.len() {
        // We use the fact that k[0] stores dydt(tn, yn).
        let error_bound = error_control.absolute_magnitude
            + error_control.relative_magnitude
                * (error_control.function_scale * y[i].abs()
                    + error_control.derivative_scale * dt * k[0][i].abs());
        max_error_ratio = f64::max(max_error_ratio, error_estimate[i] / error_bound);
    }

    (error_estimate, max_error_ratio)
}

/// Options to configure adaptive time-stepping.
pub struct AdaptiveOdeSolverOptions {
    /// Initial time for the solution.
    pub t_initial: f64,
    /// Final time for the solution.
    pub t_final: f64,
    /// Length of first attempted time step.
    pub dt_initial: f64,
    /// Parameters that determine the error bounds used to accept or reject time steps.
    pub error_control: ErrorControlOptions,
}

/// Options for computing desired error when performing adaptive time-stepping.
///
/// `rkf45_adaptive` compares its estimate to a desired error, which is computed as
/// ```
/// D = max(absolute_magnitude + relative_magnitude * (function_scale * y[i] + derivative_scale * dt * dydt[i]))
/// ```
/// where `[i]` denotes the `i`th component, and the max is taken over all components.
///
/// This form supports a variety of ways of controlling the size of the desired error
/// relative to `y` itself or to its increments.
pub struct ErrorControlOptions {
    /// Magnitude of the absolute component of desired error. Even if relative error is the primary
    /// feature of interest, this must be set to a nonzero value as a safety measure in case
    /// both y and dydt are near zero.
    pub absolute_magnitude: f64,
    /// Magnitude of the relative component of desired error.
    pub relative_magnitude: f64,
    /// Contribution of `y` to the relative component of desired error.
    pub function_scale: f64,
    /// Contribution of `y`'s increment (`dt * dydt`) to the relative component of desired error.
    pub derivative_scale: f64,
}

impl ErrorControlOptions {
    /// A simple error control option that sets the desired error to
    /// ```
    /// D = max(scale * (1 + y[i] + dt * dydt[i]))
    /// ```
    /// This has a lower bound of `scale`, but it grows proportional to `y` or `dydt` as either one
    /// becomes large.
    pub fn simple(scale: f64) -> ErrorControlOptions {
        ErrorControlOptions {
            absolute_magnitude: scale,
            relative_magnitude: scale,
            function_scale: 1.0,
            derivative_scale: 1.0,
        }
    }
}

/// Solves an initial value problem using RKF45 with adaptive time-stepping.
///
/// The method of error control and time step refinement are described in
/// Press and Teukolsky's [Adaptive Stepsize Runge-Kutta Integration, Computers in Physics 6,
/// 188 (1992)](https://doi.org/10.1063/1.4823060).
///
/// Args:
/// - `y`: Slice containing the initial value of the solution. This will be updated in-place.
/// - `dydt`: Function that evaluates the time derivative of `y`.
/// - `options`: Specifies parameters for the solver.
///
/// Returns:
/// - Vector estimating total elementwise numerical error incurred by integration.
pub fn rkf45_adaptive(
    y: &mut [f64],
    dydt: &impl Fn(f64, &[f64]) -> Vec<f64>,
    options: &AdaptiveOdeSolverOptions,
) -> Result<Vec<f64>, Error> {
    macro_rules! validate_input {
        ($x:expr) => {
            if !($x) {
                return Err(format_err!("Failed input validation: {}", stringify!($x)));
            }
        };
    }

    validate_input!(options.t_final > options.t_initial);
    validate_input!(options.dt_initial > 0.0);
    validate_input!(options.error_control.absolute_magnitude > 0.0);
    validate_input!(options.error_control.relative_magnitude >= 0.0);
    validate_input!(options.error_control.function_scale >= 0.0);
    validate_input!(options.error_control.derivative_scale >= 0.0);

    let mut t = options.t_initial;
    let mut dt = options.dt_initial;
    let mut total_error = vec![0.0; y.len()];

    let mut work = y.to_vec();
    while t < options.t_final {
        if t + dt > options.t_final {
            dt = options.t_final - t;
        }

        let (error_estimate, max_error_ratio) =
            rkf45_step(&mut work, dydt, t, dt, &options.error_control);
        use adaptive_stepping_params as params;
        if max_error_ratio < params::ERROR_RATIO_REFINING_THRESHOLD {
            // Accept the step, and update time.
            y.copy_from(&work);
            total_error.add(&error_estimate);
            t += dt;
            if max_error_ratio < params::ERROR_RATIO_COARSENING_THRESHOLD {
                // Error was particularly small, so increase step size.
                let factor = f64::min(
                    params::SAFETY_FACTOR * max_error_ratio.powf(params::COARSENING_EXPONENT),
                    params::MAX_COARSENING_FACTOR,
                );
                dt *= factor;
            }
        } else {
            // Throw out this time stemp. Revert the work array, and reduce dt.
            work.copy_from(y);
            let factor = f64::max(
                params::SAFETY_FACTOR * max_error_ratio.powf(params::REFINING_EXPONENT),
                params::MIN_REFINING_FACTOR,
            );
            dt *= factor;
        }
    }

    Ok(total_error)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;
    use test_util::{assert_gt, assert_lt};

    // rkf45_step requires ErrorControlOptions as an input. But when performing convergence
    // tests on that function, the options are meaningless. This is named to document that fact.
    static MEANINGLESS_OPTIONS: ErrorControlOptions = ErrorControlOptions {
        absolute_magnitude: 1e-8,
        relative_magnitude: 1e-8,
        function_scale: 1.0,
        derivative_scale: 1.0,
    };

    // Test rkf45_step on a first-order problem:
    //     y' = lambda*y, y(0)=1.
    // The exact solution is y(t)=exp(lambda*t).
    //
    // rkf45_step should be accurate up to O(dt^5) over a single time step (local truncation
    // error). Its error estimate should be accurate up to O(dt). We test both expectations
    // by the way the numerical errors decrease as dt is refined.
    #[test]
    fn test_first_order_problem_rkf45_step() {
        let lambda = -0.1;
        let dydt = |_t: f64, y: &[f64]| -> Vec<f64> { vec![lambda * y[0]] };
        let y_true = |t: f64| -> f64 { f64::exp(lambda * t) };

        // Record the actual error values, and the errors in estimated error values (i.e.
        // how close the estimated numerical error is to the actual numerical error) over
        // succesively-halved time steps.
        let mut actual_errors = Vec::new();
        let mut errors_in_estimated_error = Vec::new();
        for dt in &[1.0, 0.5, 0.25] {
            let mut y = [1.0];
            let estimated_error = rkf45_step(&mut y, &dydt, 0.0, *dt, &MEANINGLESS_OPTIONS).0[0];
            let actual_error = (y_true(*dt) - y[0]).abs();
            errors_in_estimated_error.push((estimated_error - actual_error).abs());
            actual_errors.push(actual_error);
        }

        // Each time dt is halved, actual_errors should shrink by roughly 1/32, and
        // errors_in_estimated_error should shrink by roughly 1/2. The factor of 1.10 allows a 10%
        // margin versus the expected convergence rate.
        assert_lt!(actual_errors[1], actual_errors[0] / 32.0 * 1.10);
        assert_lt!(actual_errors[2], actual_errors[1] / 32.0 * 1.10);
        assert_lt!(errors_in_estimated_error[1], errors_in_estimated_error[0] / 2.0 * 1.10);
        assert_lt!(errors_in_estimated_error[2], errors_in_estimated_error[1] / 2.0 * 1.10);
    }

    // Test rkf45_adaptive on the same first-order problem as above.
    //
    // This checks that we integrate to t_final with the requested degree of accuracy.
    #[test]
    fn test_first_order_problem_rkf45_adaptive() -> Result<(), Error> {
        let lambda = -0.1;
        let dydt = |_t: f64, y: &[f64]| -> Vec<f64> { vec![lambda * y[0]] };
        let y_true = |t: f64| -> f64 { f64::exp(lambda * t) };

        let options = AdaptiveOdeSolverOptions {
            t_initial: 0.0,
            t_final: 3.0,
            dt_initial: 0.1,
            error_control: ErrorControlOptions::simple(1e-6),
        };
        let mut y = [1.0];
        rkf45_adaptive(&mut y, &dydt, &options)?;
        let actual_error = (y_true(options.t_final) - y[0]).abs();

        // Error should be near the size specified, but not smaller than it.
        assert_lt!(actual_error, 1e-5);
        assert_gt!(actual_error, 1e-7);

        Ok(())
    }

    // Test rkf45_step for a second-order problem:
    //     y'' = -y; y(0)=1, y'(0)=0.
    // The exact solution is y(t)=cos(t).
    //
    // To apply the ODE solver, we translate this to the first-order system
    //     y[0]' = y[1]
    //     y[1]' = -y[0].
    // with initial conditions y[0](0)=1, y[1](0)=0.
    //
    // This test follows the same methodology as test_first_order_problem_rkf45_step;
    // see it for detailed documentation.
    #[test]
    fn test_second_order_problem_rkf45_step() {
        let dydt = |_t: f64, y: &[f64]| -> Vec<f64> { vec![y[1], -y[0]] };
        let y_true = |t: f64| -> f64 { f64::cos(t) };

        let mut actual_errors = Vec::new();
        let mut errors_in_estimated_error = Vec::new();
        for dt in &[PI / 4.0, PI / 8.0, PI / 16.0] {
            let mut y = [1.0, 0.0];
            let estimated_error = rkf45_step(&mut y, &dydt, 0.0, *dt, &MEANINGLESS_OPTIONS).0[0];
            let actual_error = (y_true(*dt) - y[0]).abs();
            errors_in_estimated_error.push((estimated_error - actual_error).abs());
            actual_errors.push(actual_error);
        }

        assert_lt!(actual_errors[1], actual_errors[0] / 32.0 * 1.10);
        assert_lt!(actual_errors[2], actual_errors[1] / 32.0 * 1.10);
        assert_lt!(errors_in_estimated_error[1], errors_in_estimated_error[0] / 2.0 * 1.10);
        assert_lt!(errors_in_estimated_error[2], errors_in_estimated_error[1] / 2.0 * 1.10);
    }

    // Test rkf45_adaptive for the same second-order problem as above.
    #[test]
    fn test_second_order_problem_rkf45_adaptive() -> Result<(), Error> {
        let dydt = |_t: f64, y: &[f64]| -> Vec<f64> { vec![y[1], -y[0]] };
        let y_true = |t: f64| -> f64 { f64::cos(t) };

        let options = AdaptiveOdeSolverOptions {
            t_initial: 0.0,
            t_final: 2.0 * PI,
            dt_initial: PI / 4.0,
            error_control: ErrorControlOptions::simple(1e-6),
        };
        let mut y = [1.0, 0.0];
        rkf45_adaptive(&mut y, &dydt, &options)?;
        let actual_error = (y_true(options.t_final) - y[0]).abs();

        // Error should be near the size specified, but not smaller than it.
        assert_lt!(actual_error, 1e-4);
        assert_gt!(actual_error, 1e-7);

        Ok(())
    }

    // Test rkf45_step for a third-order time-variant problem.
    //
    // This is a contrived problem with solution
    //     y(t) = cos(alpha*t^2).
    // As a third-order scalar equation,
    //     y''' + 4*alpha^2*t^2*y' + 12*alpha^2*t*y = 0; y(0)=1, y'(0)=0, y''(0)=0.
    // As a first-order system,
    //     y[0]' = y[1]
    //     y[1]' = y[2].
    //     y[2]' = -12*alpha^2*t*y[0] - 4*alpha^2*t^2*y[1]
    // with initial conditions y[0](0)=1, y[1](0)=0, y[2](0)=0.
    //
    // This test follows the same methodology as test_first_order_problem_rkf45_step;
    // see it for detailed documentation.
    #[test]
    fn test_third_order_problem_rkf45_step() {
        let alpha = 0.1;
        let square = |x: f64| -> f64 { x * x };
        let dydt = |t: f64, y: &[f64]| -> Vec<f64> {
            vec![y[1], y[2], -12.0 * square(alpha) * t * y[0] - 4.0 * square(alpha * t) * y[1]]
        };
        let y_true = |t: f64| -> f64 { f64::cos(alpha * square(t)) };

        let mut actual_errors = Vec::new();
        let mut errors_in_estimated_error = Vec::new();
        for dt in &[0.25, 0.125, 0.0625] {
            let mut y = [1.0, 0.0, 0.0];
            let estimated_error = rkf45_step(&mut y, &dydt, 0.0, *dt, &MEANINGLESS_OPTIONS).0[0];
            let actual_error = (y_true(*dt) - y[0]).abs();
            errors_in_estimated_error.push((estimated_error - actual_error).abs());
            actual_errors.push(actual_error);
        }

        assert_lt!(actual_errors[1], actual_errors[0] / 32.0 * 1.10);
        assert_lt!(actual_errors[2], actual_errors[1] / 32.0 * 1.10);
        assert_lt!(errors_in_estimated_error[1], errors_in_estimated_error[0] / 2.0 * 1.10);
        assert_lt!(errors_in_estimated_error[2], errors_in_estimated_error[1] / 2.0 * 1.10);
    }
    // This tests rkf45_adaptive for the same third-order problem as above.
    #[test]
    fn test_third_order_problem_rkf45_adaptive() -> Result<(), Error> {
        let alpha = 0.1;
        let square = |x: f64| -> f64 { x * x };
        let dydt = |t: f64, y: &[f64]| -> Vec<f64> {
            vec![y[1], y[2], -12.0 * square(alpha) * t * y[0] - 4.0 * square(alpha * t) * y[1]]
        };
        let y_true = |t: f64| -> f64 { f64::cos(alpha * square(t)) };

        let options = AdaptiveOdeSolverOptions {
            t_initial: 0.0,
            t_final: 4.0,
            dt_initial: 0.25,
            error_control: ErrorControlOptions::simple(1e-6),
        };
        let mut y = [1.0, 0.0, 0.0];
        rkf45_adaptive(&mut y, &dydt, &options)?;
        let actual_error = (y_true(options.t_final) - y[0]).abs();

        // Error should be near the size specified, but not smaller than it.
        assert_lt!(actual_error, 1e-4);
        assert_gt!(actual_error, 1e-7);

        Ok(())
    }

    // Test rkf45_adaptive on a problem with multiple time scales.
    //
    // This is the canonical use case for adaptive time-stepping. We have a problem with two time
    // scales, one that is much faster than the other. Our t_final is 1, and we suggest a single
    // step of length 1, which would be plenty for the slow time scale. However, the fast time
    // scale requires much smaller time steps, and the integrator must detect this.
    #[test]
    fn test_multiple_time_scales() -> Result<(), Error> {
        let lambda1 = 10.0;
        let lambda2 = 0.001;
        let dydt = |_t: f64, y: &[f64]| -> Vec<f64> { vec![-lambda1 * y[0], -lambda2 * y[1]] };
        let y_true = |t: f64| -> Vec<f64> { vec![f64::exp(-lambda1 * t), f64::exp(-lambda2 * t)] };

        let options = AdaptiveOdeSolverOptions {
            t_initial: 0.0,
            t_final: 1.0,
            dt_initial: 1.0,
            error_control: ErrorControlOptions::simple(1e-6),
        };
        let mut y = [1.0, 1.0];
        rkf45_adaptive(&mut y, &dydt, &options)?;
        let mut actual_error = y_true(options.t_final);
        actual_error.iter_mut().zip(y.iter()).for_each(|(p, q)| *p = (*p - *q).abs());

        // Our requested error scale bounds the error in both components above. The second component
        // has a much smaller error than our requested tolerance because it varies so slowly relative
        // to the first, which governs the time step selection. Hence, we don't test actual_error[1]
        // for a lower bound.
        assert_lt!(actual_error[0], 1e-5);
        assert_lt!(actual_error[1], 1e-5);
        assert_gt!(actual_error[0], 1e-7);

        Ok(())
    }

    // Test that rkf45_advance returns an error when passed invalid options.
    #[test]
    fn test_error_checks() {
        let dydt = |_t: f64, _y: &[f64]| -> Vec<f64> { vec![0.0] };

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 2.0, // Greater than t_final
                t_final: 1.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions::simple(1e-8),
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: -0.1, // Negative
                error_control: ErrorControlOptions {
                    absolute_magnitude: 1e-8,
                    relative_magnitude: 1e-8,
                    function_scale: 1.0,
                    derivative_scale: 1.0,
                }
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions {
                    absolute_magnitude: -1e-8, // Negative
                    relative_magnitude: 1e-8,
                    function_scale: 1.0,
                    derivative_scale: 1.0,
                }
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions {
                    absolute_magnitude: 1e-8,
                    relative_magnitude: -1e-8, // Negative
                    function_scale: 1.0,
                    derivative_scale: 1.0,
                }
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions {
                    absolute_magnitude: 1e-8,
                    relative_magnitude: 1e-8,
                    function_scale: -1.0, // Negative
                    derivative_scale: 1.0,
                }
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions {
                    absolute_magnitude: 1e-8,
                    relative_magnitude: 1e-8,
                    function_scale: 1.0,
                    derivative_scale: -1.0, // Negative
                }
            }
        )
        .is_err());

        assert!(rkf45_adaptive(
            &mut [1.0],
            &dydt,
            &AdaptiveOdeSolverOptions {
                t_initial: 1.0,
                t_final: 2.0,
                dt_initial: 0.1,
                error_control: ErrorControlOptions {
                    absolute_magnitude: 0.0, // Must be positive
                    relative_magnitude: 1e-8,
                    function_scale: 1.0,
                    derivative_scale: 1.0,
                }
            }
        )
        .is_err());
    }
}