Merge pull request #156 from hugohp/multiroot

Adding root-finding with derivatives to multiroot

gsl-sys/src/auto.rs

@@ -22147,7 +22147,32 @@ extern "C" {
 }
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]
-pub struct gsl_multiroot_function_fdf_struct;
+pub struct gsl_multiroot_function_fdf_struct {
+    pub f: ::std::option::Option<
+        unsafe extern "C" fn(
+            x: *const gsl_vector,
+            params: *mut ::std::os::raw::c_void,
+            f: *mut gsl_vector,
+        ) -> ::std::os::raw::c_int,
+    >,
+    pub df: ::std::option::Option<
+        unsafe extern "C" fn(
+            x: *const gsl_vector,
+            params: *mut ::std::os::raw::c_void,
+            J: *mut gsl_matrix,
+        ) -> ::std::os::raw::c_int,
+    >,
+    pub fdf: ::std::option::Option<
+        unsafe extern "C" fn(
+            x: *const gsl_vector,
+            params: *mut ::std::os::raw::c_void,
+            f: *mut gsl_vector,
+            J: *mut gsl_matrix,
+        ) -> ::std::os::raw::c_int,
+    >,
+    pub n: usize,
+    pub params: *mut ::std::os::raw::c_void,
+}
 pub type gsl_multiroot_function_fdf = gsl_multiroot_function_fdf_struct;
 #[repr(C)]
 #[derive(Debug, Copy, Clone)]

src/types/multiroot.rs

@@ -57,7 +57,7 @@ The algorithms estimate the matrix J or J^{-1} by approximate methods.
 !*/
 
 use crate::ffi::FFI;
-use crate::{Value, VectorF64, View};
+use crate::{MatrixF64, Value, VectorF64, View};
 use sys::libc::{c_int, c_void};
 
 ffi_wrapper!(
@@ -71,7 +71,7 @@ ffi_wrapper!(
 );
 
 impl MultiRootFSolverType {
-    ///This is a version of the Hybrid algorithm which replaces calls to the Jacobian function by
+    /// This is a version of the Hybrid algorithm which replaces calls to the Jacobian function by
     /// its finite difference approximation. The finite difference approximation is computed
     /// using `gsl_multiroots_fdjac()` with a relative step size of `GSL_SQRT_DBL_EPSILON`.
     /// Note that this step size will not be suitable for all problems.
@@ -234,6 +234,247 @@ impl<'a> MultiRootFSolver<'a> {
     }
 }
 
+ffi_wrapper!(
+    MultiRootFdfSolverType,
+    *const sys::gsl_multiroot_fdfsolver_type
+);
+
+impl MultiRootFdfSolverType {
+    /// This is a modified version of Powell’s Hybrid method as implemented in the HYBRJ algorithm in MINPACK.
+    /// Minpack was written by Jorge J. Moré, Burton S. Garbow and Kenneth E. Hillstrom.
+    /// The Hybrid algorithm retains the fast convergence of Newton’s method but will also reduce the
+    /// residual when Newton’s method is unreliable.
+    ///
+    /// The algorithm uses a generalized trust region to keep each step under control.
+    /// In order to be accepted a proposed new position x' must satisfy the condition |D (x' - x)| < \delta,
+    /// where D is a diagonal scaling matrix and \delta is the size of the trust region.
+    /// The components of D are computed internally, using the column norms of the Jacobian to estimate the
+    /// sensitivity of the residual to each component of x. This improves the behavior of the algorithm for badly scaled functions.
+    ///
+    /// On each iteration the algorithm first determines the standard Newton step by solving the system
+    /// J dx = - f. If this step falls inside the trust region it is used as a trial step in the
+    /// next stage. If not, the algorithm uses the linear combination of the Newton and gradient directions
+    /// which is predicted to minimize the norm of the function while staying inside the trust region,
+    /// dx = - \alpha J^{-1} f(x) - \beta \nabla |f(x)|^2
+    /// This combination of Newton and gradient directions is referred to as a dogleg step.
+    ///
+    /// The proposed step is now tested by evaluating the function at the resulting point, x'.
+    /// If the step reduces the norm of the function sufficiently then it is accepted and size of
+    /// the trust region is increased. If the proposed step fails to improve the solution then the
+    /// size of the trust region is decreased and another trial step is computed.
+    ///
+    /// The speed of the algorithm is increased by computing the changes to the Jacobian approximately,
+    /// using a rank-1 update. If two successive attempts fail to reduce the residual then the full
+    /// Jacobian is recomputed. The algorithm also monitors the progress of the solution and returns an error if several steps fail to make any improvement,
+    /// `crate::Value::NoProgress` the iteration is not making any progress, preventing the algorithm from continuing.
+    /// `crate::Value::NoProgressJacobian` re-evaluations of the Jacobian indicate that the iteration is not making any progress, preventing the algorithm from continuing.
+    #[doc(alias = "gsl_multiroot_fdfsolver_hybridsj")]
+    pub fn hybridsj() -> Self {
+        ffi_wrap!(gsl_multiroot_fdfsolver_hybridsj)
+    }
+
+    /// This algorithm is an unscaled version of `hybridsj`. The steps are controlled by a spherical trust region
+    /// |x' - x| < \delta, instead of a generalized region. This can be useful if the generalized region estimated
+    /// by `hybridsj` is inappropriate.
+    #[doc(alias = "gsl_multiroot_fdfsolver_hybridj")]
+    pub fn hybridj() -> Self {
+        ffi_wrap!(gsl_multiroot_fdfsolver_hybridj)
+    }
+
+    /// Newton’s Method is the standard root-polishing algorithm. The algorithm begins with
+    /// an initial guess for the location of the solution. On each iteration a linear approximation to
+    /// the function F is used to estimate the step which will zero all the components of the residual.
+    /// The iteration is defined by the following sequence, x \to x' = x - J^{-1} f(x)
+    /// where the Jacobian matrix J is computed from the derivative functions provided by f.
+    /// The step dx is obtained by solving the linear system, J dx = - f(x)
+    /// using LU decomposition. If the Jacobian matrix is singular, an error code of
+    /// `crate::Value::Domain` is returned.
+    #[doc(alias = "gsl_multiroot_fdfsolver_newton")]
+    pub fn newton() -> Self {
+        ffi_wrap!(gsl_multiroot_fdfsolver_newton)
+    }
+
+    /// This is a modified version of Newton’s method which attempts to improve global convergence
+    /// by requiring every step to reduce the Euclidean norm of the residual, |f(x)|.
+    /// If the Newton step leads to an increase in the norm then a reduced step of relative size,
+    /// t = (\sqrt{1 + 6 r} - 1) / (3 r)
+    /// is proposed, with r being the ratio of norms |f(x')|^2/|f(x)|^2.
+    /// This procedure is repeated until a suitable step size is found.
+    #[doc(alias = "gsl_multiroot_fdfsolver_gnewton")]
+    pub fn gnewton() -> Self {
+        ffi_wrap!(gsl_multiroot_fdfsolver_gnewton)
+    }
+}
+
+pub struct MultiRootFdfSolverFunction<'a> {
+    pub f: Box<dyn Fn(&VectorF64, &mut VectorF64) -> Value + 'a>,
+    pub df: Box<dyn Fn(&VectorF64, &mut MatrixF64) -> Value + 'a>,
+    pub fdf: Box<dyn Fn(&VectorF64, &mut VectorF64, &mut MatrixF64) -> Value + 'a>,
+    pub n: usize,
+    intern: sys::gsl_multiroot_function_fdf,
+}
+
+impl<'a> MultiRootFdfSolverFunction<'a> {
+    #[doc(alias = "gsl_multiroot_function_fdf")]
+    pub fn new<
+        F: Fn(&VectorF64, &mut VectorF64) -> Value + 'a,
+        DF: Fn(&VectorF64, &mut MatrixF64) -> Value + 'a,
+        FDF: Fn(&VectorF64, &mut VectorF64, &mut MatrixF64) -> Value + 'a,
+    >(
+        f: F,
+        df: DF,
+        fdf: FDF,
+        n: usize,
+    ) -> MultiRootFdfSolverFunction<'a> {
+        unsafe extern "C" fn inner_f(
+            x: *const sys::gsl_vector,
+            params: *mut c_void,
+            f: *mut sys::gsl_vector,
+        ) -> i32 {
+            let t = &*(params as *mut MultiRootFdfSolverFunction);
+            let i_f = &t.f;
+            i_f(
+                &VectorF64::soft_wrap(x as *const _ as *mut _),
+                &mut VectorF64::soft_wrap(f as *const _ as *mut _),
+            )
+            .into()
+        }
+
+        unsafe extern "C" fn inner_df(
+            x: *const sys::gsl_vector,
+            params: *mut c_void,
+            J: *mut sys::gsl_matrix,
+        ) -> i32 {
+            let t = &*(params as *mut MultiRootFdfSolverFunction);
+            let i_df = &t.df;
+            i_df(
+                &VectorF64::soft_wrap(x as *const _ as *mut _),
+                &mut MatrixF64::soft_wrap(J as *const _ as *mut _),
+            )
+            .into()
+        }
+
+        unsafe extern "C" fn inner_fdf(
+            x: *const sys::gsl_vector,
+            params: *mut c_void,
+            f: *mut sys::gsl_vector,
+            J: *mut sys::gsl_matrix,
+        ) -> i32 {
+            let t = &*(params as *mut MultiRootFdfSolverFunction);
+            let i_fdf = &t.fdf;
+            i_fdf(
+                &VectorF64::soft_wrap(x as *const _ as *mut _),
+                &mut VectorF64::soft_wrap(f as *const _ as *mut _),
+                &mut MatrixF64::soft_wrap(J as *const _ as *mut _),
+            )
+            .into()
+        }
+
+        MultiRootFdfSolverFunction {
+            f: Box::new(f),
+            df: Box::new(df),
+            fdf: Box::new(fdf),
+            n,
+            intern: sys::gsl_multiroot_function_fdf {
+                f: Some(inner_f),
+                df: Some(inner_df),
+                fdf: Some(inner_fdf),
+                n,
+                params: std::ptr::null_mut(),
+            },
+        }
+    }
+
+    #[allow(clippy::wrong_self_convention)]
+    fn to_raw(&mut self) -> *mut sys::gsl_multiroot_function_fdf {
+        self.intern.n = self.n;
+        self.intern.params = self as *mut MultiRootFdfSolverFunction as *mut c_void;
+        &mut self.intern
+    }
+}
+
+ffi_wrapper!(
+    MultiRootFdfSolver,
+    *mut sys::gsl_multiroot_fdfsolver,
+    gsl_multiroot_fdfsolver_free
+);
+
+impl MultiRootFdfSolver {
+    /// This function returns a pointer to a newly allocated instance of a derivative solver of type T
+    /// for a system of `n` dimensions.
+    #[doc(alias = "gsl_multiroot_fdfsolver_alloc")]
+    pub fn new(t: MultiRootFdfSolverType, n: usize) -> Option<MultiRootFdfSolver> {
+        let ptr = unsafe { sys::gsl_multiroot_fdfsolver_alloc(t.unwrap_shared(), n) };
+
+        if ptr.is_null() {
+            None
+        } else {
+            Some(Self::wrap(ptr))
+        }
+    }
+
+    /// These functions set, or reset, an existing solver to use the functions `f`, and the initial guess `x.
+    #[doc(alias = "gsl_multiroot_fdfsolver_set")]
+    pub fn set(&mut self, f: &mut MultiRootFdfSolverFunction, x: &VectorF64) -> Result<(), Value> {
+        let ret = unsafe {
+            sys::gsl_multiroot_fdfsolver_set(self.unwrap_unique(), f.to_raw(), x.unwrap_shared())
+        };
+        result_handler!(ret, ())
+    }
+
+    /// Return the name of the solver.
+    #[doc(alias = "gsl_multiroot_fdfsolver_name")]
+    pub fn name(&self) -> Option<String> {
+        let n = unsafe { sys::gsl_multiroot_fdfsolver_name(self.unwrap_shared()) };
+        if n.is_null() {
+            return None;
+        }
+        let mut len = 0;
+        loop {
+            if unsafe { *n.offset(len) } == 0 {
+                break;
+            }
+            len += 1;
+        }
+        let slice = unsafe { std::slice::from_raw_parts(n as _, len as _) };
+        std::str::from_utf8(slice).ok().map(|x| x.to_owned())
+    }
+
+    /// Perform a single iteration of the solver. If the iteration encounters an
+    /// unexpected problem then an error code will be returned,
+    ///
+    /// * `crate::Value::BadFunc` the iteration encountered a singular point where the function or its derivative evaluated to Inf or NaN.
+    ///
+    /// * `crate::Value::NoProgress` the iteration is not making any progress,
+    /// preventing the algorithm from continuing.
+    ///
+    /// The solver maintains a current best estimate of the root and its function value at all times.
+    /// This information can be accessed with `root`, `f`, and `dx` functions.
+    #[doc(alias = "gsl_multiroot_fdfsolver_iterate")]
+    pub fn iterate(&mut self) -> Result<(), Value> {
+        let ret = unsafe { sys::gsl_multiroot_fdfsolver_iterate(self.unwrap_unique()) };
+        result_handler!(ret, ())
+    }
+
+    /// Returns the current estimate of the root for the solver.
+    #[doc(alias = "gsl_multiroot_fdfsolver_root")]
+    pub fn root(&self) -> View<'_, VectorF64> {
+        unsafe { View::new(sys::gsl_multiroot_fdfsolver_root(self.unwrap_shared())) }
+    }
+
+    /// Returns the last step taken by the solver.
+    #[doc(alias = "gsl_multiroot_fdfsolver_dx")]
+    pub fn dx(&self) -> View<'_, VectorF64> {
+        unsafe { View::new(sys::gsl_multiroot_fdfsolver_dx(self.unwrap_shared())) }
+    }
+
+    /// Returns the function value f(x) at the current estimate of the root for the solver.
+    #[doc(alias = "gsl_multiroot_fdfsolver_f")]
+    pub fn f(&self) -> View<'_, VectorF64> {
+        unsafe { View::new(sys::gsl_multiroot_fdfsolver_f(self.unwrap_shared())) }
+    }
+}
+
 #[cfg(any(test, doctest))]
 mod tests {
     /// This doc block will be used to ensure that the closure can't be set everywhere!
@@ -277,11 +518,27 @@ mod tests {
     use super::*;
     use crate::multiroot::test_residual;
 
+    const RPARAMS: (f64, f64) = (1.0, 10.0);
+
     /// checking a test function
     /// must return a success criteria (or failure)
     fn rosenbrock_f(x: &VectorF64, f: &mut VectorF64) -> Value {
-        f.set(0, 1.0 - x.get(0));
-        f.set(1, x.get(0) - x.get(1).powf(2.0));
+        f.set(0, RPARAMS.0 * (1.0 - x.get(0)));
+        f.set(1, RPARAMS.1 * (x.get(1) - x.get(0).powf(2.0)));
+        Value::Success
+    }
+
+    fn rosenbrock_df(x: &VectorF64, J: &mut MatrixF64) -> Value {
+        J.set(0, 0, -RPARAMS.0);
+        J.set(0, 1, 0f64);
+        J.set(1, 0, -2.0 * RPARAMS.1 * x.get(0));
+        J.set(1, 1, RPARAMS.1);
+        Value::Success
+    }
+
+    fn rosenbrock_fdf(x: &VectorF64, f: &mut VectorF64, J: &mut MatrixF64) -> Value {
+        rosenbrock_f(x, f);
+        rosenbrock_df(x, J);
         Value::Success
     }
 
@@ -298,6 +555,19 @@ mod tests {
         )
     }
 
+    fn print_fdf_state(solver: &mut MultiRootFdfSolver, iteration: usize) {
+        let f = solver.f();
+        let x = solver.root();
+        println!(
+            "iter: {}, f = [{:+.2e}, {:+.2e}], x = [{:+.5}, {:+.5}]",
+            iteration,
+            f.get(0),
+            f.get(1),
+            x.get(0),
+            x.get(1)
+        )
+    }
+
     #[test]
     fn test_multiroot_fsolver() {
         // setup workspace
@@ -338,4 +608,49 @@ mod tests {
         }
         assert!(matches!(status, crate::Value::Success))
     }
+
+    #[test]
+    fn test_multiroot_fdf_fsolver() {
+        // setup workspace
+        let mut multi_root = MultiRootFdfSolver::new(MultiRootFdfSolverType::gnewton(), 2).unwrap();
+        let array_size: usize = 2;
+        let guess_value = VectorF64::from_slice(&[-10.0, -5.0]).unwrap();
+        let mut fs = MultiRootFdfSolverFunction::new(
+            rosenbrock_f,
+            rosenbrock_df,
+            rosenbrock_fdf,
+            array_size,
+        );
+        multi_root.set(&mut fs, &guess_value).unwrap();
+
+        // iteration counters
+        let max_iter: usize = 100;
+        let mut iter = 0;
+
+        // convergence checks
+        let mut status = crate::Value::Continue;
+        let epsabs = 1e-6;
+
+        print_fdf_state(&mut multi_root, 0);
+
+        while matches!(status, crate::Value::Continue) && iter < max_iter {
+            // iterate solver
+            multi_root.iterate().unwrap();
+
+            // print current iteration
+            print_fdf_state(&mut multi_root, iter);
+
+            // test for convergence
+            let f_value = multi_root.f();
+            status = test_residual(&f_value, epsabs);
+
+            // check if iteration succeeded
+            if matches!(status, crate::Value::Success) {
+                println!("Converged");
+            }
+
+            iter += 1;
+        }
+        assert!(matches!(status, crate::Value::Success))
+    }
 }