Added output function of optimum.

src/first_order.rs

@@ -171,7 +171,7 @@ where
 /// let gradf = |x: &Array1<f32>| array![2. * (x[0] + 2.) * (2. * x[0].powi(2) - x[0] - 1.)];
 /// let x0 = &array![-1.];
 ///
-/// let x_star = descent(
+/// let (x_star, f_star) = descent(
 ///     f,
 ///     gradf,
 ///     &x0,
@@ -181,7 +181,7 @@ where
 /// ).unwrap();
 /// assert!((-2. - x_star[0]).abs() < 1e-10);
 ///
-/// let x_star = descent(
+/// let (x_star, f_star) = descent(
 ///     f,
 ///     gradf,
 ///     &x0,
@@ -192,11 +192,11 @@ where
 /// assert!((-2. - x_star[0]).abs() < 1e-10);
 ///
 /// let x0 = &array![-0.5];
-/// let x_star = descent(f, gradf, &x0, Default::default(), 1e-3, Some(10)).unwrap();
+/// let (x_star, f_star) = descent(f, gradf, &x0, Default::default(), 1e-3, Some(10)).unwrap();
 /// assert!((-0.5 - x_star[0]).abs() < 1e-10);
 ///
 /// let x0 = &array![0.];
-/// let x_star = descent(f, gradf, &x0, Default::default(), 1e-3, Some(10)).unwrap();
+/// let (x_star, f_star) = descent(f, gradf, &x0, Default::default(), 1e-3, Some(10)).unwrap();
 /// assert!((1. - x_star[0]).abs() < 1e-10);
 ///
 /// // It also takes multivariate objective functions
@@ -209,9 +209,9 @@ where
 ///     ]
 /// };
 /// let x = array![1f32, -0.5f32];
-/// let opt = descent(f, gradf, &x, Default::default(), 1e-3, Some(10000)).unwrap();
-/// assert!((opt[0] - 1.).abs() <= 1e-2);
-/// assert!((opt[1] - 1.).abs() <= 1e-2);
+/// let (x_star, f_star) = descent(f, gradf, &x, Default::default(), 1e-3, Some(10000)).unwrap();
+/// assert!((x_star[0] - 1.).abs() <= 1e-2);
+/// assert!((x_star[1] - 1.).abs() <= 1e-2);
 /// ```
 pub fn descent<X, F, G>(
     f: F,
@@ -220,7 +220,7 @@ pub fn descent<X, F, G>(
     params: DescentParameter<X::Elem>,
     gtol: X::Elem,
     maxiter: Option<usize>,
-) -> Result<X, TuutalError<X>>
+) -> Result<(X, X::Elem), TuutalError<X>>
 where
     X: Vector + Clone + VecDot<Output = X::Elem>,
     for<'a> &'a X: Add<X, Output = X> + Mul<&'a X, Output = X> + Mul<X, Output = X>,

src/first_order/macros.rs

@@ -104,6 +104,7 @@ macro_rules! impl_optimizer_descent {
         {
             type Iterate = X;
             type Intermediate = $step;
+            type ObjectiveOutput = X::Elem;
             fn nb_iter(&self) -> usize {
                 self.counter.iter
             }
@@ -113,6 +114,11 @@ macro_rules! impl_optimizer_descent {
             fn intermediate(&self) -> Self::Intermediate {
                 self.sigma.clone()
             }
+            fn objective_output(&mut self) -> Self::ObjectiveOutput {
+                let fx = self.func(&self.x);
+                self.counter.fcalls += 1;
+                fx
+            }
         }
     };
 }

src/first_order/steepest_descent/unit_test.rs

@@ -20,7 +20,7 @@ mod tests {
         let armijo = DescentParameter::new_armijo(0.01, 0.5);
         let (f, gradf) = rosenbrock_2d();
         let x = Array1::from_iter([1f32, -0.5f32]);
-        let opt = descent(f, gradf, &x, armijo, 1e-4, Some(10000)).unwrap();
+        let (opt, _f_star) = descent(f, gradf, &x, armijo, 1e-4, Some(10000)).unwrap();
         let expected = Array1::from_iter([1., 1.]);
         assert!(l2_diff(&opt, &expected) < 1e-3);
     }
@@ -30,7 +30,7 @@ mod tests {
         let powolf = DescentParameter::new_powell_wolfe(0.0001, 0.9);
         let (f, gradf) = rosenbrock_2d();
         let x = Array1::from_iter([1f32, -0.5f32]);
-        let opt = descent(f, gradf, &x, powolf, 1e-4, Some(10000)).unwrap();
+        let (opt, _f_star) = descent(f, gradf, &x, powolf, 1e-4, Some(10000)).unwrap();
         let expected = Array1::from_iter([1., 1.]);
         assert!(l2_diff(&opt, &expected) < 1e-3);
     }
@@ -92,7 +92,7 @@ mod tests {
             ])
         };
         let x = Array1::from_iter([10f32, -15., -100.]);
-        let opt = descent(f, gradf, &x, powolf, 1e-4, Some(10000)).unwrap();
+        let (opt, _f_star) = descent(f, gradf, &x, powolf, 1e-4, Some(10000)).unwrap();
         let expected = Array1::from_iter([1., 1., 1.]);
         assert!(l2_diff(&opt, &expected) < 1e-3);
     }

src/traits.rs

@@ -101,6 +101,7 @@ where
 pub trait Optimizer: core::iter::Iterator {
     type Iterate;
     type Intermediate;
+    type ObjectiveOutput;
     /// Number of iterations done so far.
     fn nb_iter(&self) -> usize;
     /// Current iterate.
@@ -109,7 +110,7 @@ pub trait Optimizer: core::iter::Iterator {
     fn optimize(
         &mut self,
         maxiter: Option<usize>,
-    ) -> Result<Self::Iterate, TuutalError<Self::Iterate>> {
+    ) -> Result<(Self::Iterate, Self::ObjectiveOutput), TuutalError<Self::Iterate>> {
         let maxiter = maxiter.unwrap_or(1000);
         while let Some(_) = self.next() {
             if self.nb_iter() > maxiter {
@@ -119,10 +120,12 @@ pub trait Optimizer: core::iter::Iterator {
                 });
             }
         }
-        Ok(self.iterate())
+        Ok((self.iterate(), self.objective_output()))
     }
-    /// Gives intermediate values during the algorithm like step size.
+    /// Gives current intermediate values during the algorithm like step size.
     fn intermediate(&self) -> Self::Intermediate;
+    /// Gives objective function output for current iterate.
+    fn objective_output(&mut self) -> Self::ObjectiveOutput;
 }
 
 /// Implements the notion of upper and lower bounds

src/zero_order/nelder_mead.rs

@@ -17,15 +17,15 @@ type SimplexParameterResult<A> = Result<(A, A, A, A), TuutalError<Array1<A>>>;
 /// // Example from python scipy.optimize.minimize_scalar
 /// let f = |x: &Array1<f32>| (x[0] - 2.) * x[0] * (x[0] + 2.).powi(2);
 /// let x0 = &array![-1.];
-/// let x_star =
+/// let (x_star, f_star) =
 ///     nelder_mead::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, true, None)
 ///     .unwrap();
 /// assert!((-2. - x_star[0]).abs() <= 2e-4);
 ///
 /// let f =
 ///     |arr: &Array1<f32>| 100. * (arr[1] - arr[0].powi(2)).powi(2) + (1. - arr[0]).powi(2);
 /// let x0 = array![1., -0.5];
-/// let x_star =
+/// let (x_star, f_star) =
 ///     nelder_mead::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, true, None)
 ///     .unwrap();
 /// assert!((1. - x_star[0]).abs() <= 1e-3);
@@ -41,7 +41,7 @@ pub fn nelder_mead<A, B, F>(
     fatol: A,
     adaptive: bool,
     bounds: Option<B>,
-) -> Result<Array1<A>, TuutalError<Array1<A>>>
+) -> Result<(Array1<A>, A), TuutalError<Array1<A>>>
 where
     A: Scalar<Array1<A>>,
     B: Bound<A>,
@@ -401,7 +401,7 @@ impl<F, A> NelderMeadIterates<F, A> {
     }
 
     /// Computes the objective function value for a given input vector.
-    pub fn obj(&self, x: &Array1<A>) -> A
+    pub fn func(&self, x: &Array1<A>) -> A
     where
         F: Fn(&Array1<A>) -> A,
     {
@@ -519,7 +519,6 @@ where
     type Item = Array1<A>;
     fn next(&mut self) -> Option<Self::Item> {
         if self.convergence {
-            self.iter += 1;
             return None; // TODO
         } else if self.stop() {
             self.convergence = true;
@@ -533,7 +532,7 @@ where
         let one = A::cast_from_f32(1.);
         let xbar = self.centroid().unwrap(); // Is it allways safe to .unwrap()?
         let xr = self.affine(&xbar, self.rho, one);
-        let fxr = self.obj(&xr);
+        let fxr = self.func(&xr);
         self.fcalls += 1;
         let mut doshrink = false;
         let last = self.sim.nrows() - 1;
@@ -543,7 +542,7 @@ where
                 return None; // TODO
             }
             let xe = self.affine(&xbar, self.rho, self.chi);
-            let fxe = self.obj(&xe);
+            let fxe = self.func(&xe);
             self.fcalls += 1;
             if fxe < fxr {
                 self.sim.row_mut(last).assign(&xe);
@@ -561,7 +560,7 @@ where
                 return None; // TODO
             }
             let xc = self.affine(&xbar, self.psi, self.rho);
-            let fxc = self.obj(&xc);
+            let fxc = self.func(&xc);
             self.fcalls += 1;
             if fxc <= fxr {
                 self.sim.row_mut(last).assign(&xc);
@@ -575,7 +574,7 @@ where
         } else {
             // Perform an inside contraction
             let xcc = self.affine(&xbar, self.psi, -one);
-            let fxcc = self.obj(&xcc);
+            let fxcc = self.func(&xcc);
             self.fcalls += 1;
             if fxcc < self.fsim[last] {
                 self.sim.row_mut(last).assign(&xcc);
@@ -607,11 +606,17 @@ where
 {
     type Iterate = Array1<A>;
     type Intermediate = ();
+    type ObjectiveOutput = A;
     fn nb_iter(&self) -> usize {
         self.iter
     }
     fn iterate(&self) -> Array1<A> {
         self.sim.row(0).to_owned()
     }
     fn intermediate(&self) {}
+    fn objective_output(&mut self) -> Self::ObjectiveOutput {
+        let fx = self.func(&self.iterate());
+        self.fcalls += 1;
+        fx
+    }
 }

src/zero_order/nelder_mead/unit_test.rs

@@ -125,10 +125,10 @@ mod test {
     fn test_nelder_mead() {
         let f = |arr: &Array1<f32>| arr.dot(arr);
         let x0 = Array1::from_iter([-5., -5.]);
-        let x_star =
+        let (x_star, f_star) =
             nelder_mead::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, true, None)
                 .unwrap();
         assert!(l2_diff(&x_star, &Array1::from_iter([0., 0.])) < 2e-3);
-        assert!(f(&x_star) < 1e-5);
+        assert!(f_star < 1e-5);
     }
 }

src/zero_order/powell.rs

@@ -18,15 +18,15 @@ use super::{default_nb_iter, scalar::BrentOptResult};
 /// // Example from python scipy.optimize.minimize_scalar
 /// let f = |x: &Array1<f32>| (x[0] - 2.) * x[0] * (x[0] + 2.).powi(2);
 /// let x0 = &array![-1.];
-/// let x_star =
+/// let (x_star, f_star) =
 ///     powell::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, None)
 ///     .unwrap();
 /// assert!((x_star[0] - 1.280776).abs() <= 1e-4);
 ///
 /// let f =
 ///     |arr: &Array1<f32>| 100. * (arr[1] - arr[0].powi(2)).powi(2) + (1. - arr[0]).powi(2);
 /// let x0 = array![1., -0.5];
-/// let x_star =
+/// let (x_star, f_star) =
 ///     powell::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, None)
 ///     .unwrap();
 /// assert!((x_star[0] - 1.).abs() <= 1e-5);
@@ -41,7 +41,7 @@ pub fn powell<A, B, F>(
     xtol: A,
     ftol: A,
     bounds: Option<B>,
-) -> Result<Array1<A>, TuutalError<Array1<A>>>
+) -> Result<(Array1<A>, A), TuutalError<Array1<A>>>
 where
     A: Scalar<Array1<A>> + core::fmt::Debug,
     B: Bound<A>,
@@ -272,7 +272,7 @@ impl<F, A> PowellIterates<F, A> {
         })
     }
 
-    pub(crate) fn obj(&self, x: &Array1<A>) -> A
+    pub(crate) fn func(&self, x: &Array1<A>) -> A
     where
         F: Fn(&Array1<A>) -> A,
     {
@@ -355,7 +355,7 @@ where
             self.iter += 1;
             return None; // TO change
         }
-        let fx2 = self.obj(&x2);
+        let fx2 = self.func(&x2);
         self.fcalls += 1;
         if fx > fx2 {
             let mut t = two * (fx + fx2 - two * self.fval);
@@ -401,11 +401,17 @@ where
 {
     type Iterate = Array1<A>;
     type Intermediate = ();
+    type ObjectiveOutput = A;
     fn nb_iter(&self) -> usize {
         self.iter
     }
     fn iterate(&self) -> Array1<A> {
         self.x.clone()
     }
     fn intermediate(&self) {}
+    fn objective_output(&mut self) -> Self::ObjectiveOutput {
+        let fx = self.func(&self.iterate());
+        self.fcalls += 1;
+        fx
+    }
 }

src/zero_order/powell/unit_test.rs

@@ -49,9 +49,9 @@ mod tests {
     fn test_powell() {
         let f = |arr: &Array1<f32>| arr.dot(arr);
         let x0 = Array1::from_iter([-5., -5.]);
-        let x_star =
+        let (x_star, f_star) =
             powell::<_, (f32, f32), _>(f, &x0, None, Some(100), None, 1e-5, 1e-5, None).unwrap();
         assert!(l2_diff(&x_star, &Array1::from_iter([0., 0.])) < 1e-6);
-        assert!(f(&x_star) < 1e-6);
+        assert!(f_star < 1e-6);
     }
 }