feat!: Add Function and Optimizer traits for general optimization

Although the primary goal of Gomez - solving non-linear systems of
equations - remains still the same, I believe that providing means for
general numerical optimization may benefit someone for the features that
Gomez has: first-class support for bounds, simple, low-level, iterative
API and various utilities. Moreover, two out of three algorithms
implemented so far are optimization algorithms, not solvers anyway.

As solving equations is the primary goal, the information about general
optimization is rather omitted from the documentation and that is on
purpose.

Author: pnevyk

benches/solvers.rs

@@ -33,7 +33,7 @@ where
     let mut fx = x.clone_owned();
     let mut iter = 0;
     loop {
-        if solver.next(&f, &dom, &mut x, &mut fx).is_err() {
+        if solver.next(f, dom, &mut x, &mut fx).is_err() {
             return false;
         }
 
@@ -67,10 +67,7 @@ fn rosenbrock1(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -95,10 +92,7 @@ fn rosenbrock2(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -123,10 +117,7 @@ fn rosenbrock_scaled(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -147,10 +138,7 @@ fn rosenbrock_large(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -175,10 +163,7 @@ fn powell(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -199,10 +184,7 @@ fn bullard_biegler(c: &mut Criterion) {
             assert!(solve(
                 &f,
                 &dom,
-                GslSolverWrapper::new(GslFunctionWrapper::new(
-                    f.clone(),
-                    GslVec::from(x.as_slice())
-                )),
+                GslSolverWrapper::new(GslFunctionWrapper::new(f, GslVec::from(x.as_slice()))),
                 x.clone_owned()
             ))
         })
@@ -249,7 +231,7 @@ where
             na::U1::name(),
         );
 
-        match self.f.apply(&x, &mut fx) {
+        match self.f.eval(&x, &mut fx) {
             Ok(_) => GslStatus::ok(),
             Err(_) => GslStatus::err(GslError::BadFunc),
         }

examples/rosenbrock.rs

@@ -9,15 +9,21 @@ struct Rosenbrock {
     b: f64,
 }
 
-impl System for Rosenbrock {
+impl Problem for Rosenbrock {
     type Scalar = f64;
     type Dim = na::U2;
 
-    fn apply<Sx, Sfx>(
+    fn dim(&self) -> Self::Dim {
+        na::U2::name()
+    }
+}
+
+impl System for Rosenbrock {
+    fn eval<Sx, Sfx>(
         &self,
         x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
         fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
-    ) -> Result<(), SystemError>
+    ) -> Result<(), Error>
     where
         Sx: na::storage::Storage<Self::Scalar, Self::Dim>,
         Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
@@ -27,10 +33,6 @@ impl System for Rosenbrock {
 
         Ok(())
     }
-
-    fn dim(&self) -> Self::Dim {
-        na::U2::name()
-    }
 }
 
 fn main() -> Result<(), String> {

src/analysis/initial.rs

@@ -16,7 +16,7 @@ use nalgebra::{
 use thiserror::Error;
 
 use crate::{
-    core::{Domain, System, SystemError},
+    core::{Domain, Error, Problem, System},
     derivatives::{Jacobian, JacobianError, EPSILON_SQRT},
 };
 
@@ -25,14 +25,14 @@ use crate::{
 pub enum InitialGuessAnalysisError {
     /// Error that occurred when evaluating the system.
     #[error("{0}")]
-    System(#[from] SystemError),
+    System(#[from] Error),
     /// Error that occurred when computing the Jacobian matrix.
     #[error("{0}")]
     Jacobian(#[from] JacobianError),
 }
 
 /// Initial guesses analyzer. See [module](self) documentation for more details.
-pub struct InitialGuessAnalysis<F: System> {
+pub struct InitialGuessAnalysis<F: Problem> {
     nonlinear: Vec<usize>,
     ty: PhantomData<F>,
 }
@@ -59,7 +59,7 @@ where
         let scale = OVector::from_iterator_generic(f.dim(), U1::name(), scale_iter);
 
         // Compute F'(x) in the initial point.
-        f.apply(x, fx)?;
+        f.eval(x, fx)?;
         let jac1 = Jacobian::new(f, x, &scale, fx)?;
 
         // Compute Newton step.
@@ -74,7 +74,7 @@ where
         *x += p;
 
         // Compute F'(x) after one Newton step.
-        f.apply(x, fx)?;
+        f.eval(x, fx)?;
         let jac2 = Jacobian::new(f, x, &scale, fx)?;
 
         // Linear variables have no effect on the Jacobian matrix. They can be

src/core.rs

@@ -1,17 +1,23 @@
 //! Core abstractions and types for Gomez.
 //!
-//! *Users* are mainly interested in implementing the [`System`] trait,
-//! optionally specifying the [domain](Domain).
+//! *Users* are mainly interested in implementing the [`System`] and [`Problem`]
+//! traits, optionally specifying the [domain](Domain).
 //!
 //! Algorithms *developers* are interested in implementing the [`Solver`] trait
-//! and using extension traits [`SystemExt`] and [`VectorDomainExt`] as well as
-//! tools in [derivatives](crate::derivatives) or
-//! [population](crate::population) modules.
+//! (or possibly the [`Optimizer`] trait too) and using extension traits (e.g.,
+//! [`VectorDomainExt`]) as well as tools in [derivatives](crate::derivatives)
+//! or [population](crate::population) modules.
 
+mod base;
 mod domain;
+mod function;
+mod optimizer;
 mod solver;
 mod system;
 
+pub use base::*;
 pub use domain::*;
+pub use function::*;
+pub use optimizer::*;
 pub use solver::*;
 pub use system::*;

src/core/base.rs

@@ -0,0 +1,42 @@
+use nalgebra::{Dim, RealField};
+use thiserror::Error;
+
+use super::domain::Domain;
+
+/// The base trait for [`System`](super::system::System) and
+/// [`Function`](super::function::Function).
+pub trait Problem {
+    /// Type of the scalar, usually f32 or f64.
+    type Scalar: RealField;
+
+    /// Dimension of the system. Can be fixed
+    /// ([`Const`](nalgebra::base::dimension::Const)) or dynamic
+    /// ([`Dynamic`](nalgebra::base::dimension::Dynamic)).
+    type Dim: Dim;
+
+    /// Return the actual dimension of the system. This is needed for dynamic
+    /// systems.
+    fn dim(&self) -> Self::Dim;
+
+    /// Get the domain (bound constraints) of the system. If not overridden, the
+    /// system is unconstrained.
+    fn domain(&self) -> Domain<Self::Scalar> {
+        Domain::with_dim(self.dim().value())
+    }
+}
+
+/// Error encountered while applying variables to the function.
+#[derive(Debug, Error)]
+pub enum Error {
+    /// The number of variables does not match the dimensionality
+    /// ([`Problem::dim`]) of the problem.
+    #[error("invalid dimensionality")]
+    InvalidDimensionality,
+    /// An invalid value (NaN, positive or negative infinity) of a residual or
+    /// the function value occurred.
+    #[error("invalid value encountered")]
+    InvalidValue,
+    /// A custom error specific to the system or function.
+    #[error("{0}")]
+    Custom(Box<dyn std::error::Error>),
+}

src/core/function.rs

@@ -0,0 +1,109 @@
+use nalgebra::{
+    allocator::Allocator, storage::Storage, storage::StorageMut, DefaultAllocator, Vector,
+};
+use num_traits::Zero;
+
+use super::{
+    base::{Error, Problem},
+    system::System,
+};
+
+/// The trait for defining functions.
+///
+/// ## Defining a function
+///
+/// A function is any type that implements [`Function`] and [`Problem`] traits.
+/// There are two required associated types (scalar type and dimension type) and
+/// two required methods: [`apply`](Function::apply) and [`dim`](Problem::dim).
+///
+/// ```rust
+/// use gomez::nalgebra as na;
+/// use gomez::prelude::*;
+/// use na::{Dim, DimName};
+///
+/// // A problem is represented by a type.
+/// struct Rosenbrock {
+///     a: f64,
+///     b: f64,
+/// }
+///
+/// impl Problem for Rosenbrock {
+///     // The numeric type. Usually f64 or f32.
+///     type Scalar = f64;
+///     // The dimension of the problem. Can be either statically known or dynamic.
+///     type Dim = na::U2;
+///
+///     // Return the actual dimension of the system.
+///     fn dim(&self) -> Self::Dim {
+///         na::U2::name()
+///     }
+/// }
+///
+/// impl Function for Rosenbrock {
+///     // Apply trial values of variables to the function.
+///     fn apply<Sx>(
+///         &self,
+///         x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
+///     ) -> Result<Self::Scalar, Error>
+///     where
+///         Sx: na::storage::Storage<Self::Scalar, Self::Dim>,
+///     {
+///         // Compute the function value.
+///         Ok((self.a - x[0]).powi(2) + self.b * (x[1] - x[0].powi(2)).powi(2))
+///     }
+/// }
+/// ```
+pub trait Function: Problem {
+    /// Calculate the function value given values of the variables.
+    fn apply<Sx>(&self, x: &Vector<Self::Scalar, Self::Dim, Sx>) -> Result<Self::Scalar, Error>
+    where
+        Sx: Storage<Self::Scalar, Self::Dim>;
+
+    /// Calculate the norm of residuals of the system given values of the
+    /// variable for cases when the function is actually a system of equations.
+    ///
+    /// The optimizers should prefer calling this function because the
+    /// implementation for systems reuse `fx` for calculating the residuals and
+    /// do not make an unnecessary allocation for it.
+    fn apply_eval<Sx, Sfx>(
+        &self,
+        x: &Vector<Self::Scalar, Self::Dim, Sx>,
+        fx: &mut Vector<Self::Scalar, Self::Dim, Sfx>,
+    ) -> Result<Self::Scalar, Error>
+    where
+        Sx: Storage<Self::Scalar, Self::Dim>,
+        Sfx: StorageMut<Self::Scalar, Self::Dim>,
+    {
+        let norm = self.apply(x)?;
+        fx.fill(Self::Scalar::zero());
+        fx[0] = norm;
+        Ok(norm)
+    }
+}
+
+impl<F: System> Function for F
+where
+    DefaultAllocator: Allocator<F::Scalar, F::Dim>,
+{
+    fn apply<Sx>(&self, x: &Vector<Self::Scalar, Self::Dim, Sx>) -> Result<Self::Scalar, Error>
+    where
+        Sx: Storage<Self::Scalar, Self::Dim>,
+    {
+        let mut fx = x.clone_owned();
+        self.apply_eval(x, &mut fx)
+    }
+
+    fn apply_eval<Sx, Sfx>(
+        &self,
+        x: &Vector<Self::Scalar, Self::Dim, Sx>,
+        fx: &mut Vector<Self::Scalar, Self::Dim, Sfx>,
+    ) -> Result<Self::Scalar, Error>
+    where
+        Sx: Storage<Self::Scalar, Self::Dim>,
+        Sfx: StorageMut<Self::Scalar, Self::Dim>,
+    {
+        self.eval(x, fx)?;
+        let norm = fx.norm();
+        Ok(norm)
+    }
+}