feat!: introduce own RealField trait

src/analysis/initial.rs

@@ -14,8 +14,8 @@ use nalgebra::{
 };
 
 use crate::{
-    core::{Domain, Problem, System},
-    derivatives::{Jacobian, EPSILON_SQRT},
+    core::{Domain, Problem, RealField, System},
+    derivatives::Jacobian,
 };
 
 /// Initial guesses analyzer. See [module](self) documentation for more details.
@@ -75,7 +75,7 @@ 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();

src/core/base.rs

@@ -1,7 +1,18 @@
-use nalgebra::RealField;
-
 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 {
@@ -12,3 +23,13 @@ pub trait Problem {
     /// system is unconstrained.
     fn domain(&self) -> Domain<Self::Field>;
 }
+
+impl RealField for f32 {
+    const EPSILON_SQRT: Self = 0.00034526698;
+    const EPSILON_CBRT: Self = 0.0049215667;
+}
+
+impl RealField for f64 {
+    const EPSILON_SQRT: Self = 0.000000014901161193847656;
+    const EPSILON_CBRT: Self = 0.0000060554544523933395;
+}

src/derivatives.rs

@@ -3,21 +3,12 @@
 use std::ops::Deref;
 
 use nalgebra::{
-    convert,
     storage::{Storage, StorageMut},
     ComplexField, DimName, Dynamic, IsContiguous, OMatrix, OVector, RealField, Vector, U1,
 };
 use num_traits::{One, Zero};
 
-use crate::core::{Function, Problem, System};
-
-/// Square root of double precision machine epsilon. This value is a standard
-/// constant for epsilons in approximating first-order derivate-based concepts.
-pub const EPSILON_SQRT: f64 = 0.000000014901161193847656;
-
-/// Cubic root of double precision machine epsilon. This value is a standard
-/// constant for epsilons in approximating second-order derivate-based concepts.
-pub const EPSILON_CBRT: f64 = 0.0000060554544523933395;
+use crate::core::{Function, Problem, RealField as _, System};
 
 /// Jacobian matrix of a system.
 #[derive(Debug)]
@@ -76,7 +67,7 @@ impl<F: System> Jacobian<F> {
         Sscale: Storage<F::Field, Dynamic>,
         Sfx: Storage<F::Field, Dynamic>,
     {
-        let eps: F::Field = convert(EPSILON_SQRT);
+        let eps = F::Field::EPSILON_SQRT;
 
         for (j, mut col) in self.jac.column_iter_mut().enumerate() {
             let xj = x[j];
@@ -173,7 +164,7 @@ impl<F: Function> Gradient<F> {
         Sx: StorageMut<F::Field, Dynamic> + IsContiguous,
         Sscale: Storage<F::Field, Dynamic>,
     {
-        let eps: F::Field = convert(EPSILON_SQRT);
+        let eps = F::Field::EPSILON_SQRT;
 
         for i in 0..x.nrows() {
             let xi = x[i];
@@ -265,7 +256,7 @@ impl<F: Function> Hessian<F> {
         Sx: StorageMut<F::Field, Dynamic> + IsContiguous,
         Sscale: Storage<F::Field, Dynamic>,
     {
-        let eps: F::Field = convert(EPSILON_CBRT);
+        let eps = F::Field::EPSILON_CBRT;
 
         for i in 0..x.nrows() {
             let xi = x[i];

src/solver/nelder_mead.rs

@@ -31,10 +31,7 @@ use nalgebra::{
 use num_traits::{One, Zero};
 use thiserror::Error;
 
-use crate::{
-    core::{Domain, Function, Optimizer, Problem, Solver, System},
-    derivatives::EPSILON_SQRT,
-};
+use crate::core::{Domain, Function, Optimizer, Problem, RealField as _, Solver, System};
 
 /// Family of coefficients for reflection, expansion and contractions.
 #[derive(Debug, Clone, Copy)]
@@ -435,7 +432,7 @@ impl<F: Function> NelderMead<F> {
             // otherwise an error reduction was achieved. This criterion is
             // taken from "Less is more: Simplified Nelder-Mead method for large
             // unconstrained optimization".
-            let eps: F::Field = convert(EPSILON_SQRT);
+            let eps = F::Field::EPSILON_SQRT;
 
             let worst = errors[sort_perm[n]];
             let best = errors[sort_perm[0]];

src/solver/trust_region.rs

@@ -37,8 +37,8 @@ use num_traits::{One, Zero};
 use thiserror::Error;
 
 use crate::{
-    core::{Domain, Function, Optimizer, Problem, Solver, System},
-    derivatives::{Gradient, Hessian, Jacobian, EPSILON_SQRT},
+    core::{Domain, Function, Optimizer, Problem, RealField as _, Solver, System},
+    derivatives::{Gradient, Hessian, Jacobian},
 };
 
 /// Specification for initial value of trust region size.
@@ -102,7 +102,7 @@ impl<F: Problem> TrustRegionOptions<F> {
 impl<F: Problem> Default for TrustRegionOptions<F> {
     fn default() -> Self {
         Self {
-            delta_min: convert(EPSILON_SQRT),
+            delta_min: F::Field::EPSILON_SQRT,
             delta_max: convert(1e9),
             delta_init: DeltaInit::Estimated,
             mu_min: convert(1e-10),