API improvements

examples/rosenbrock.rs

@@ -1,6 +1,6 @@
 use gomez::nalgebra as na;
-use gomez::prelude::*;
 use gomez::solver::TrustRegion;
+use gomez::{Domain, Problem, Solver, System};
 use na::{Dynamic, IsContiguous};
 
 // https://en.wikipedia.org/wiki/Rosenbrock_function
@@ -10,27 +10,24 @@ struct Rosenbrock {
 }
 
 impl Problem for Rosenbrock {
-    type Scalar = f64;
+    type Field = f64;
 
-    fn domain(&self) -> Domain<Self::Scalar> {
+    fn domain(&self) -> Domain<Self::Field> {
         Domain::unconstrained(2)
     }
 }
 
 impl System for Rosenbrock {
     fn eval<Sx, Sfx>(
         &self,
-        x: &na::Vector<Self::Scalar, Dynamic, Sx>,
-        fx: &mut na::Vector<Self::Scalar, Dynamic, Sfx>,
-    ) -> Result<(), ProblemError>
-    where
-        Sx: na::storage::Storage<Self::Scalar, Dynamic> + IsContiguous,
-        Sfx: na::storage::StorageMut<Self::Scalar, Dynamic>,
+        x: &na::Vector<Self::Field, Dynamic, Sx>,
+        fx: &mut na::Vector<Self::Field, Dynamic, Sfx>,
+    ) where
+        Sx: na::storage::Storage<Self::Field, Dynamic> + IsContiguous,
+        Sfx: na::storage::StorageMut<Self::Field, Dynamic>,
     {
         fx[0] = (self.a - x[0]).powi(2);
         fx[1] = self.b * (x[1] - x[0].powi(2)).powi(2);
-
-        Ok(())
     }
 }
 
@@ -46,7 +43,7 @@ fn main() -> Result<(), String> {
 
     for i in 1..=100 {
         solver
-            .next(&f, &dom, &mut x, &mut fx)
+            .solve_next(&f, &dom, &mut x, &mut fx)
             .map_err(|err| format!("{}", err))?;
 
         println!(

gomez-bench/benches/main.rs

@@ -5,9 +5,9 @@ fn main() {
 use gomez::{
     nalgebra as na,
     nalgebra::Dynamic,
-    prelude::*,
     solver::{NelderMead, TrustRegion},
     testing::*,
+    Domain, Problem, Solver,
 };
 use gsl_wrapper::{
     function::Function as GslFunction,
@@ -168,7 +168,7 @@ mod bullard_biegler {
 fn bench_solve<F, S, GF, GS>(bencher: divan::Bencher, with_system: GF, with_solver: GS)
 where
     GF: Fn() -> (F, usize),
-    GS: Fn(&F, &Domain<F::Scalar>, &na::OVector<F::Scalar, Dynamic>) -> S,
+    GS: Fn(&F, &Domain<F::Field>, &na::OVector<F::Field, Dynamic>) -> S,
     F: TestSystem,
     S: Solver<F>,
 {
@@ -184,7 +184,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.solve_next(&f, &dom, &mut x, &mut fx).is_err() {
                     panic!("solver error");
                 }
 
@@ -203,8 +203,8 @@ where
 
 fn with_trust_region<F>(
     f: &F,
-    dom: &Domain<F::Scalar>,
-    _: &na::OVector<F::Scalar, Dynamic>,
+    dom: &Domain<F::Field>,
+    _: &na::OVector<F::Field, Dynamic>,
 ) -> TrustRegion<F>
 where
     F: Problem,
@@ -214,8 +214,8 @@ where
 
 fn with_nelder_mead<F>(
     f: &F,
-    dom: &Domain<F::Scalar>,
-    _: &na::OVector<F::Scalar, Dynamic>,
+    dom: &Domain<F::Field>,
+    _: &na::OVector<F::Field, Dynamic>,
 ) -> NelderMead<F>
 where
     F: Problem,
@@ -225,11 +225,11 @@ where
 
 fn with_gsl_hybrids<F>(
     f: &F,
-    _: &Domain<F::Scalar>,
-    x: &na::OVector<F::Scalar, Dynamic>,
+    _: &Domain<F::Field>,
+    x: &na::OVector<F::Field, Dynamic>,
 ) -> GslSolverWrapper<GslFunctionWrapper<F>>
 where
-    F: TestSystem<Scalar = f64> + Clone,
+    F: TestSystem<Field = f64> + Clone,
 {
     GslSolverWrapper::new(GslFunctionWrapper::new(
         f.clone(),
@@ -248,7 +248,7 @@ impl<F> GslFunctionWrapper<F> {
     }
 }
 
-impl<F: TestSystem<Scalar = f64>> GslFunction for GslFunctionWrapper<F> {
+impl<F: TestSystem<Field = f64>> GslFunction for GslFunctionWrapper<F> {
     fn eval(&self, x: &GslVec, f: &mut GslVec) -> GslStatus {
         use na::DimName;
         let dim = Dynamic::new(x.len());
@@ -264,10 +264,8 @@ impl<F: TestSystem<Scalar = f64>> GslFunction for GslFunctionWrapper<F> {
             na::U1::name(),
         );
 
-        match self.f.eval(&x, &mut fx) {
-            Ok(_) => GslStatus::ok(),
-            Err(_) => GslStatus::err(GslError::BadFunc),
-        }
+        self.f.eval(&x, &mut fx);
+        GslStatus::ok()
     }
 
     fn init(&self) -> GslVec {
@@ -287,21 +285,21 @@ impl<F: GslFunction> GslSolverWrapper<F> {
     }
 }
 
-impl<F: TestSystem<Scalar = f64>> Solver<F> for GslSolverWrapper<GslFunctionWrapper<F>> {
+impl<F: TestSystem<Field = f64>> Solver<F> for GslSolverWrapper<GslFunctionWrapper<F>> {
     const NAME: &'static str = "GSL hybrids";
 
     type Error = String;
 
-    fn next<Sx, Sfx>(
+    fn solve_next<Sx, Sfx>(
         &mut self,
         _f: &F,
-        _dom: &Domain<F::Scalar>,
-        x: &mut na::Vector<F::Scalar, Dynamic, Sx>,
-        fx: &mut na::Vector<F::Scalar, Dynamic, Sfx>,
+        _dom: &Domain<F::Field>,
+        x: &mut na::Vector<F::Field, Dynamic, Sx>,
+        fx: &mut na::Vector<F::Field, Dynamic, Sfx>,
     ) -> Result<(), Self::Error>
     where
-        Sx: na::storage::StorageMut<F::Scalar, Dynamic> + IsContiguous,
-        Sfx: na::storage::StorageMut<F::Scalar, Dynamic>,
+        Sx: na::storage::StorageMut<F::Field, Dynamic> + IsContiguous,
+        Sfx: na::storage::StorageMut<F::Field, Dynamic>,
     {
         let result = self.solver.step().to_result();
         x.copy_from_slice(self.solver.root());

src/analysis/initial.rs

@@ -12,24 +12,12 @@ use std::marker::PhantomData;
 use nalgebra::{
     convert, storage::StorageMut, ComplexField, DimName, Dynamic, IsContiguous, OVector, Vector, U1,
 };
-use thiserror::Error;
 
 use crate::{
-    core::{Domain, Problem, ProblemError, System},
-    derivatives::{Jacobian, JacobianError, EPSILON_SQRT},
+    core::{Domain, Problem, RealField, System},
+    derivatives::Jacobian,
 };
 
-/// Error returned from [`InitialGuessAnalysis`] solver.
-#[derive(Debug, Error)]
-pub enum InitialGuessAnalysisError {
-    /// Error that occurred when evaluating the system.
-    #[error("{0}")]
-    System(#[from] ProblemError),
-    /// 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: Problem> {
     nonlinear: Vec<usize>,
@@ -40,13 +28,13 @@ impl<F: System> InitialGuessAnalysis<F> {
     /// Analyze the system in given point.
     pub fn analyze<Sx, Sfx>(
         f: &F,
-        dom: &Domain<F::Scalar>,
-        x: &mut Vector<F::Scalar, Dynamic, Sx>,
-        fx: &mut Vector<F::Scalar, Dynamic, Sfx>,
-    ) -> Result<Self, InitialGuessAnalysisError>
+        dom: &Domain<F::Field>,
+        x: &mut Vector<F::Field, Dynamic, Sx>,
+        fx: &mut Vector<F::Field, Dynamic, Sfx>,
+    ) -> Self
     where
-        Sx: StorageMut<F::Scalar, Dynamic> + IsContiguous,
-        Sfx: StorageMut<F::Scalar, Dynamic>,
+        Sx: StorageMut<F::Field, Dynamic> + IsContiguous,
+        Sfx: StorageMut<F::Field, Dynamic>,
     {
         let dim = Dynamic::new(dom.dim());
         let scale = dom
@@ -55,8 +43,8 @@ impl<F: System> InitialGuessAnalysis<F> {
             .unwrap_or_else(|| OVector::from_element_generic(dim, U1::name(), convert(1.0)));
 
         // Compute F'(x) in the initial point.
-        f.eval(x, fx)?;
-        let jac1 = Jacobian::new(f, x, &scale, fx)?;
+        f.eval(x, fx);
+        let jac1 = Jacobian::new(f, x, &scale, fx);
 
         // Compute Newton step.
         let mut p = fx.clone_owned();
@@ -66,12 +54,12 @@ impl<F: System> InitialGuessAnalysis<F> {
         qr.solve_mut(&mut p);
 
         // Do Newton step.
-        p *= convert::<_, F::Scalar>(0.001);
+        p *= convert::<_, F::Field>(0.001);
         *x += p;
 
         // Compute F'(x) after one Newton step.
-        f.eval(x, fx)?;
-        let jac2 = Jacobian::new(f, x, &scale, 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
         // recognized by observing no change in corresponding columns (i.e.,
@@ -87,15 +75,15 @@ impl<F: System> InitialGuessAnalysis<F> {
             .filter(|(_, (c1, c2))| {
                 c1.iter()
                     .zip(c2.iter())
-                    .any(|(a, b)| (*a - *b).abs() > convert(EPSILON_SQRT))
+                    .any(|(a, b)| (*a - *b).abs() > F::Field::EPSILON_SQRT)
             })
             .map(|(col, _)| col)
             .collect();
 
-        Ok(Self {
+        Self {
             nonlinear,
             ty: PhantomData,
-        })
+        }
     }
 
     /// Returns indices of variables that have influence on the Jacobian matrix

src/core/base.rs

@@ -1,31 +1,35 @@
-use nalgebra::RealField;
-use thiserror::Error;
-
 use super::domain::Domain;
 
+/// Trait implemented by real numbers.
+pub trait RealField: nalgebra::RealField {
+    /// Square root of double precision machine epsilon. This value is a
+    /// standard constant for epsilons in approximating first-order
+    /// derivate-based concepts.
+    const EPSILON_SQRT: Self;
+
+    /// Cubic root of double precision machine epsilon. This value is a standard
+    /// constant for epsilons in approximating second-order derivate-based
+    /// concepts.
+    const EPSILON_CBRT: Self;
+}
+
 /// 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 + Copy;
+    /// Type of the field, usually f32 or f64.
+    type Field: RealField + Copy;
 
     /// Get the domain (bound constraints) of the system. If not overridden, the
     /// system is unconstrained.
-    fn domain(&self) -> Domain<Self::Scalar>;
+    fn domain(&self) -> Domain<Self::Field>;
+}
+
+impl RealField for f32 {
+    const EPSILON_SQRT: Self = 0.00034526698;
+    const EPSILON_CBRT: Self = 0.0049215667;
 }
 
-/// Error encountered while applying variables to the problem.
-#[derive(Debug, Error)]
-pub enum ProblemError {
-    /// The number of variables does not match the dimensionality of the problem
-    /// domain.
-    #[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>),
+impl RealField for f64 {
+    const EPSILON_SQRT: Self = 0.000000014901161193847656;
+    const EPSILON_CBRT: Self = 0.0000060554544523933395;
 }