Added some more optimization tests.

--HG--
extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%403927

dlib/test/CMakeLists.txt

@@ -61,6 +61,7 @@ set (tests
    metaprogramming.cpp
    multithreaded_object.cpp
    optimization.cpp
+   optimization_test_functions.cpp
    opt_qp_solver.cpp
    pipe.cpp
    pixel.cpp

dlib/test/makefile

@@ -71,6 +71,7 @@ SRC += member_function_pointer.cpp
 SRC += metaprogramming.cpp
 SRC += multithreaded_object.cpp
 SRC += optimization.cpp
+SRC += optimization_test_functions.cpp
 SRC += opt_qp_solver.cpp
 SRC += pipe.cpp
 SRC += pixel.cpp

dlib/test/optimization_test_functions.cpp

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#include "optimization_test_functions.h"
+
+/*
+
+    Most of the code in this file is converted from the set of Fortran 90 routines 
+    created by John Burkardt.
+
+    The original Fortran can be found here: http://orion.math.iastate.edu/burkardt/f_src/testopt/testopt.html
+
+*/
+
+
+namespace dlib
+{
+    namespace test_functions
+    {
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,0,1> chebyquad_residuals(const matrix<double,0,1>& x)
+        {
+            matrix<double,0,1> fvec(x.size());
+            const int n = x.size();
+            int i;
+            int j;
+            double t;
+            double t1;
+            double t2;
+            double th;
+            fvec = 0;
+
+            for (j = 1; j <= n; ++j)
+            {
+                t1 = 1.0E+00;
+                t2 = 2.0E+00 * x(j-1) - 1.0E+00;
+                t = 2.0E+00 * t2;
+                for (i = 1; i <= n; ++i)
+                {
+                    fvec(i-1) = fvec(i-1) + t2;
+                    th = t * t2 - t1;
+                    t1 = t2;
+                    t2 = th;
+                }
+            }
+
+            for (i = 1; i <= n; ++i)
+            {
+                fvec(i-1) = fvec(i-1) / (double) ( n );
+                if ( ( i%2 ) == 0 ) 
+                    fvec(i-1) = fvec(i-1) + 1.0E+00 / ( (double)i*i - 1.0E+00 );
+            }
+
+            return fvec;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        double chebyquad_residual(int i, const matrix<double,0,1>& x)
+        {
+            return chebyquad_residuals(x)(i);
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        int& chebyquad_calls() 
+        {
+            static int count = 0;
+            return count;
+        }
+
+        double chebyquad(const matrix<double,0,1>& x )
+        {
+            chebyquad_calls()++;
+            return sum(squared(chebyquad_residuals(x)));
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,0,1> chebyquad_derivative (const matrix<double,0,1>& x)
+        {
+            const int n = x.size();
+            matrix<double,0,1> fvec = chebyquad_residuals(x);
+            matrix<double,0,1> g(n);
+            int i;
+            int j;
+            double s1;
+            double s2;
+            double t;
+            double t1;
+            double t2;
+            double th;
+
+            for (j = 1; j <= n; ++j)
+            {
+                g(j-1) = 0.0E+00;
+                t1 = 1.0E+00;
+                t2 = 2.0E+00 * x(j-1) - 1.0E+00;
+                t = 2.0E+00 * t2;
+                s1 = 0.0E+00;
+                s2 = 2.0E+00;
+                for (i = 1; i <= n; ++i)
+                {
+                    g(j-1) = g(j-1) + fvec(i-1) * s2;
+                    th = 4.0E+00 * t2 + t * s2 - s1;
+                    s1 = s2;
+                    s2 = th;
+                    th = t * t2 - t1;
+                    t1 = t2;
+                    t2 = th;
+                }
+            }
+
+            g = 2.0E+00 * g / (double) ( n );
+
+            return g;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,0,1> chebyquad_start (int n)
+        {
+            int i;
+            matrix<double,0,1> x(n);
+
+            for (i = 1; i <= n; ++i)
+                x(i-1) = double ( i ) / double ( n + 1 );
+
+            return x;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,0,1> chebyquad_solution (int n)
+        {
+            matrix<double,0,1> x(n);
+
+            x = 0;
+            switch (n)
+            {
+                case 2:
+                    x = 0.2113249E+00, 0.7886751E+00;
+                    break;
+                case 4:
+                    x = 0.1026728E+00, 0.4062037E+00, 0.5937963E+00, 0.8973272E+00;
+                    break;
+                case 6:
+                    x = 0.066877E+00, 0.288741E+00, 0.366682E+00, 0.633318E+00, 0.711259E+00, 0.933123E+00;
+                    break;
+                case 8:
+                    x = 0.043153E+00, 0.193091E+00, 0.266329E+00, 0.500000E+00, 0.500000E+00, 0.733671E+00, 0.806910E+00, 0.956847E+00;
+                    break;
+                default:
+                    std::ostringstream sout;
+                    sout << "don't know chebyquad solution for n = " << n;
+                    throw dlib::error(sout.str());
+                    break;
+            }
+
+            return x;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double> chebyquad_hessian(const matrix<double,0,1>& x)
+        {
+            const int lda = x.size();
+            const int n = x.size();
+            double d1;
+            double d2;
+            matrix<double,0,1> fvec = chebyquad_residuals(x);
+            matrix<double,0,1> gvec(n);
+            matrix<double> h(lda,n);
+            int i;
+            int j;
+            int k;
+            double p1;
+            double p2;
+            double s1;
+            double s2;
+            double ss1;
+            double ss2;
+            double t;
+            double t1;
+            double t2;
+            double th;
+            double tt;
+            double tth;
+            double tt1;
+            double tt2;
+            h = 0;
+
+            d1 = 1.0E+00 / double ( n );
+            d2 = 2.0E+00 * d1;
+
+            for (j = 1; j <= n; ++j)
+            {
+
+                h(j-1,j-1) = 4.0E+00 * d1;
+                t1 = 1.0E+00;
+                t2 = 2.0E+00 * x(j-1) - 1.0E+00;
+                t = 2.0E+00 * t2;
+                s1 = 0.0E+00;
+                s2 = 2.0E+00;
+                p1 = 0.0E+00;
+                p2 = 0.0E+00;
+                gvec(0) = s2;
+
+                for (i = 2; i <= n; ++i)
+                {
+                    th = 4.0E+00 * t2 + t * s2 - s1;
+                    s1 = s2;
+                    s2 = th;
+                    th = t * t2 - t1;
+                    t1 = t2;
+                    t2 = th;
+                    th = 8.0E+00 * s1 + t * p2 - p1;
+                    p1 = p2;
+                    p2 = th;
+                    gvec(i-1) = s2;
+                    h(j-1,j-1) = h(j-1,j-1) + fvec(i-1) * th + d1 * s2*s2;
+                }
+
+                h(j-1,j-1) = d2 * h(j-1,j-1);
+
+                for (k = 1; k <= j-1; ++k)
+                {
+
+                    h(j-1,k-1) = 0.0;
+                    tt1 = 1.0E+00;
+                    tt2 = 2.0E+00 * x(k-1) - 1.0E+00;
+                    tt = 2.0E+00 * tt2;
+                    ss1 = 0.0E+00;
+                    ss2 = 2.0E+00;
+
+                    for (i = 1; i <= n; ++i)
+                    {
+                        h(j-1,k-1) = h(j-1,k-1) + ss2 * gvec(i-1);
+                        tth = 4.0E+00 * tt2 + tt * ss2 - ss1;
+                        ss1 = ss2;
+                        ss2 = tth;
+                        tth = tt * tt2 - tt1;
+                        tt1 = tt2;
+                        tt2 = tth;
+                    }
+
+                    h(j-1,k-1) = d2 * d1 * h(j-1,k-1);
+
+                }
+
+            }
+
+            h = make_symmetric(h);
+            return h;
+        }
+
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+
+        double brown_residual (int i, const matrix<double,4,1>& x)
+        /*!
+            requires
+                - 1 <= i <= 20
+            ensures
+                - returns the ith brown residual
+        !*/
+        {
+            double c;
+            double f;
+            double f1;
+            double f2;
+
+            f = 0.0E+00;
+
+
+            c = double ( i ) / 5.0E+00;
+            f1 = x(0) + c * x(1) - std::exp ( c );
+            f2 = x(2) + std::sin ( c ) * x(3) - std::cos ( c );
+
+            f = f1*f1 + f2*f2; 
+
+            return f;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        double brown ( const matrix<double,4,1>& x)
+        {
+            double f;
+            int i;
+
+            f = 0;
+
+            for (i = 1; i <= 20; ++i)
+            {
+                f += std::pow(brown_residual(i, x), 2);
+            }
+
+            return f;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,4,1> brown_derivative ( const matrix<double,4,1>& x)
+        {
+            double c;
+            double df1dx1;
+            double df1dx2;
+            double df2dx3;
+            double df2dx4;
+            double f1;
+            double f2;
+            matrix<double,4,1> g;
+            int i;
+
+            g = 0;
+
+            for (i = 1; i <= 20; ++i)
+            {
+
+                c = double ( i ) / 5.0E+00;
+
+                f1 = x(0) + c * x(1) - std::exp ( c );
+                f2 = x(2) + std::sin ( c ) * x(3) - std::cos ( c );
+
+                df1dx1 = 1.0E+00;
+                df1dx2 = c;
+                df2dx3 = 1.0E+00;
+                df2dx4 = std::sin ( c );
+
+                using std::pow;
+                g(0) = g(0) + 4.0E+00 * ( pow(f1,3) * df1dx1 + f1 * pow(f2,2) * df1dx1 );
+                g(1) = g(1) + 4.0E+00 * ( pow(f1,3) * df1dx2 + f1 * pow(f2,2) * df1dx2 );
+                g(2) = g(2) + 4.0E+00 * ( pow(f1,2) * f2 * df2dx3 + pow(f2,3) * df2dx3 );
+                g(3) = g(3) + 4.0E+00 * ( pow(f1,2) * f2 * df2dx4 + pow(f2,3) * df2dx4 );
+
+            }
+
+            return g;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,4,4> brown_hessian ( const matrix<double,4,1>& x)
+        {
+            double c;
+            double df1dx1;
+            double df1dx2;
+            double df2dx3;
+            double df2dx4;
+            double f1;
+            double f2;
+            matrix<double,4,4> h;
+            int i;
+
+            h = 0;
+
+            for (i = 1; i <= 20; ++i)
+            {
+
+                c = double ( i ) / 5.0E+00;
+
+                f1 = x(0) + c * x(1) - std::exp ( c );
+                f2 = x(2) + std::sin ( c ) * x(3) - std::cos ( c );
+
+                df1dx1 = 1.0E+00;
+                df1dx2 = c;
+                df2dx3 = 1.0E+00;
+                df2dx4 = std::sin ( c );
+
+                using std::pow;
+                h(0,0) = h(0,0) + 12.0E+00 * pow(f1,2) * df1dx1 * df1dx1 + 4.0E+00 * pow(f2,2) * df1dx1 * df1dx1;
+                h(0,1) = h(0,1) + 12.0E+00 * pow(f1,2) * df1dx1 * df1dx2 + 4.0E+00 * pow(f2,2) * df1dx1 * df1dx2;
+                h(0,2) = h(0,2) + 8.0E+00 * f1 * f2 * df1dx1 * df2dx3;
+                h(0,3) = h(0,3) + 8.0E+00 * f1 * f2 * df1dx1 * df2dx4;
+
+                h(1,0) = h(1,0) + 12.0E+00 * pow(f1,2) * df1dx2 * df1dx1 + 4.0E+00 * pow(f2,2) * df1dx2 * df1dx1;
+                h(1,1) = h(1,1) + 12.0E+00 * pow(f1,2) * df1dx2 * df1dx2 + 4.0E+00 * pow(f2,2) * df1dx2 * df1dx1;
+                h(1,2) = h(1,2) + 8.0E+00 * f1 * f2 * df1dx2 * df2dx3;
+                h(1,3) = h(1,3) + 8.0E+00 * f1 * f2 * df1dx2 * df2dx4;
+
+                h(2,0) = h(2,0) + 8.0E+00 * f1 * f2 * df2dx3 * df1dx1;
+                h(2,1) = h(2,1) + 8.0E+00 * f1 * f2 * df2dx3 * df1dx2;
+                h(2,2) = h(2,2) + 4.0E+00 * pow(f1,2) * df2dx3 * df2dx3 + 12.0E+00 * pow(f2,2) * df2dx3 * df2dx3;
+                h(2,3) = h(2,3) + 4.0E+00 * pow(f1,2) * df2dx4 * df2dx3 + 12.0E+00 * pow(f2,2) * df2dx3 * df2dx4;
+
+                h(3,0) = h(3,0) + 8.0E+00 * f1 * f2 * df2dx4 * df1dx1;
+                h(3,1) = h(3,1) + 8.0E+00 * f1 * f2 * df2dx4 * df1dx2;
+                h(3,2) = h(3,2) + 4.0E+00 * pow(f1,2) * df2dx3 * df2dx4 + 12.0E+00 * pow(f2,2) * df2dx4 * df2dx3;
+                h(3,3) = h(3,3) + 4.0E+00 * pow(f1,2) * df2dx4 * df2dx4 + 12.0E+00 * pow(f2,2) * df2dx4 * df2dx4;
+
+            }
+
+            return make_symmetric(h);
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,4,1> brown_start ()
+        {
+            matrix<double,4,1> x;
+            x = 25.0E+00, 5.0E+00, -5.0E+00, -1.0E+00;
+            return x;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,4,1> brown_solution ()
+        {
+            matrix<double,4,1> x;
+            // solution from original documentation.
+            //x = -11.5844E+00, 13.1999E+00, -0.406200E+00, 0.240998E+00;
+            x = -11.594439905669450042, 13.203630051593080452, -0.40343948856573402795, 0.23677877338218666914;
+            return x;
+        }
+
+    // ----------------------------------------------------------------------------------------
+
+    }
+}
+
+

dlib/test/optimization_test_functions.h

+// Copyright (C) 2010  Davis E. King (davis@dlib.net)
+// License: Boost Software License   See LICENSE.txt for the full license.
+#ifndef DLIB_OPTIMIZATION_TEST_FUNCTiONS_H___
+#define DLIB_OPTIMIZATION_TEST_FUNCTiONS_H___
+
+#include <dlib/matrix.h>
+#include <sstream>
+#include <cmath>
+
+/*
+
+    Most of the code in this file is converted from the set of Fortran 90 routines 
+    created by John Burkardt.
+
+    The original Fortran can be found here: http://orion.math.iastate.edu/burkardt/f_src/testopt/testopt.html
+
+*/
+
+
+namespace dlib
+{
+    namespace test_functions
+    {
+
+    // ----------------------------------------------------------------------------------------
+
+        matrix<double,0,1> chebyquad_residuals(const matrix<double,0,1>& x);
+
+        double chebyquad_residual(int i, const matrix<double,0,1>& x);
+
+        int& chebyquad_calls();
+
+        double chebyquad(const matrix<double,0,1>& x );
+
+        matrix<double,0,1> chebyquad_derivative (const matrix<double,0,1>& x);
+
+        matrix<double,0,1> chebyquad_start (int n);
+
+        matrix<double,0,1> chebyquad_solution (int n);
+
+        matrix<double> chebyquad_hessian(const matrix<double,0,1>& x);
+
+    // ----------------------------------------------------------------------------------------
+
+        class chebyquad_function_model 
+        {
+        public:
+
+            // Define the type used to represent column vectors
+            typedef matrix<double,0,1> column_vector;
+            // Define the type used to represent the hessian matrix
+            typedef matrix<double> general_matrix;
+
+            double operator() ( 
+                const column_vector& x
+            ) const
+            {
+                return chebyquad(x);
+            }
+
+            void get_derivative_and_hessian (
+                const column_vector& x,
+                column_vector& d,
+                general_matrix& h
+            ) const
+            {
+                d = chebyquad_derivative(x);
+                h = chebyquad_hessian(x);
+            }
+        };
+
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+
+        double brown_residual (int i, const matrix<double,4,1>& x);
+        /*!
+            requires
+                - 1 <= i <= 20
+            ensures
+                - returns the ith brown residual
+        !*/
+
+        double brown ( const matrix<double,4,1>& x);
+
+        matrix<double,4,1> brown_derivative ( const matrix<double,4,1>& x);
+
+        matrix<double,4,4> brown_hessian ( const matrix<double,4,1>& x);
+
+        matrix<double,4,1> brown_start ();
+
+        matrix<double,4,1> brown_solution ();
+
+        class brown_function_model 
+        {
+        public:
+
+            // Define the type used to represent column vectors
+            typedef matrix<double,4,1> column_vector;
+            // Define the type used to represent the hessian matrix
+            typedef matrix<double> general_matrix;
+
+            double operator() ( 
+                const column_vector& x
+            ) const
+            {
+                return brown(x);
+            }
+
+            void get_derivative_and_hessian (
+                const column_vector& x,
+                column_vector& d,
+                general_matrix& h
+            ) const
+            {
+                d = brown_derivative(x);
+                h = brown_hessian(x);
+            }
+        };
+
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+
+        template <typename T>
+        matrix<T,2,1> rosen_big_start()
+        {
+            matrix<T,2,1> x;
+            x = -1.2, -1;
+            return x;
+        }
+
+    // This is a variation on the Rosenbrock test function but with large residuals.  The
+    // minimum is at 1, 1 and the objective value is 1.
+        template <typename T>
+        T rosen_big_residual (int i, const matrix<T,2,1>& m)
+        {
+            using std::pow;
+            const T x = m(0); 
+            const T y = m(1);
+
+            if (i == 1)
+            {
+                return 100*pow(y - x*x,2)+1.0;
+            }
+            else 
+            {
+                return pow(1 - x,2) + 1.0;
+            }
+        }
+
+        template <typename T>
+        T rosen_big ( const matrix<T,2,1>& m)
+        {
+            using std::pow;
+            return 0.5*(pow(rosen_big_residual(1,m),2) + pow(rosen_big_residual(2,m)));
+        }
+
+        template <typename T>
+        matrix<T,2,1> rosen_big_solution ()
+        {
+            matrix<T,2,1> x;
+            // solution from original documentation.
+            x = 1,1;
+            return x;
+        }
+
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+    // ----------------------------------------------------------------------------------------
+
+        template <typename T>
+        matrix<T,2,1> rosen_start()
+        {
+            matrix<T,2,1> x;
+            x = -1.2, -1;
+            return x;
+        }
+
+        template <typename T>
+        T rosen ( const matrix<T,2,1>& m)
+        {
+            const T x = m(0); 
+            const T y = m(1);
+
+            using std::pow;
+            // compute Rosenbrock's function and return the result
+            return 100.0*pow(y - x*x,2) + pow(1 - x,2);
+        }
+
+        template <typename T>
+        T rosen_residual (int i, const matrix<T,2,1>& m)
+        {
+            const T x = m(0); 
+            const T y = m(1);
+
+
+            if (i == 1)
+            {
+                return 10*(y - x*x);
+            }
+            else
+            {
+                return 1 - x;
+            }
+        }
+
+        template <typename T>
+        const matrix<T,2,1> rosen_derivative ( const matrix<T,2,1>& m)
+        {
+            const T x = m(0);
+            const T y = m(1);
+
+            // make us a column vector of length 2
+            matrix<T,2,1> res(2);
+
+            // now compute the gradient vector
+            res(0) = -400*x*(y-x*x) - 2*(1-x); // derivative of rosen() with respect to x
+            res(1) = 200*(y-x*x);              // derivative of rosen() with respect to y
+            return res;
+        }
+
+        template <typename T>
+        const matrix<T,2,2> rosen_hessian ( const matrix<T,2,1>& m)
+        {
+            const T x = m(0);
+            const T y = m(1);
+
+            // make us a column vector of length 2
+            matrix<T,2,2> res;
+
+            // now compute the gradient vector
+            res(0,0) = -400*y + 3*400*x*x + 2; 
+            res(1,1) = 200;              
+
+            res(0,1) = -400*x;              
+            res(1,0) = -400*x;              
+            return res;
+        }
+
+        template <typename T>
+        matrix<T,2,1> rosen_solution ()
+        {
+            matrix<T,2,1> x;
+            // solution from original documentation.
+            x = 1,1;
+            return x;
+        }
+
+    // ------------------------------------------------------------------------------------
+
+        template <typename T>
+        struct rosen_function_model
+        {
+            typedef matrix<T,2,1> column_vector;
+            typedef matrix<T,2,2> general_matrix;
+
+            T operator() ( column_vector x) const
+            {
+                return static_cast<T>(rosen(x));
+            }
+
+            void get_derivative_and_hessian (
+                const column_vector& x,
+                column_vector& d,
+                general_matrix& h
+            ) const 
+            {
+                d = rosen_derivative(x);
+                h = rosen_hessian(x);
+            }
+
+        };
+
+    // ----------------------------------------------------------------------------------------
+
+    }
+}
+
+#endif // DLIB_OPTIMIZATION_TEST_FUNCTiONS_H___
+
+
+

dlib/test/trust_region.cpp

@@ -3,6 +3,7 @@
 
 
 #include <dlib/optimization.h>
+#include "optimization_test_functions.h"
 #include <sstream>
 #include <string>
 #include <cstdlib>
@@ -19,79 +20,10 @@ namespace
     using namespace test;
     using namespace dlib;
     using namespace std;
+    using namespace dlib::test_functions;
 
     logger dlog("test.trust_region");
 
-// ----------------------------------------------------------------------------------------
-
-    template <typename T>
-    T rosen ( const matrix<T,2,1>& m)
-    {
-        const T x = m(0); 
-        const T y = m(1);
-
-        // compute Rosenbrock's function and return the result
-        return 100.0*pow(y - x*x,2) + pow(1 - x,2);
-    }
-
-    template <typename T>
-    const matrix<T,2,1> rosen_derivative ( const matrix<T,2,1>& m)
-    {
-        const T x = m(0);
-        const T y = m(1);
-
-        // make us a column vector of length 2
-        matrix<T,2,1> res(2);
-
-        // now compute the gradient vector
-        res(0) = -400*x*(y-x*x) - 2*(1-x); // derivative of rosen() with respect to x
-        res(1) = 200*(y-x*x);              // derivative of rosen() with respect to y
-        return res;
-    }
-
-    template <typename T>
-    const matrix<T,2,2> rosen_hessian ( const matrix<T,2,1>& m)
-    {
-        const T x = m(0);
-        const T y = m(1);
-
-        // make us a column vector of length 2
-        matrix<T,2,2> res;
-
-        // now compute the gradient vector
-        res(0,0) = -400*y + 3*400*x*x + 2; 
-        res(1,1) = 200;              
-
-        res(0,1) = -400*x;              
-        res(1,0) = -400*x;              
-        return res;
-    }
-
-// ----------------------------------------------------------------------------------------
-
-    template <typename T>
-    struct rosen_model
-    {
-        typedef matrix<T,2,1> column_vector;
-        typedef matrix<T,2,2> general_matrix;
-
-        T operator() ( column_vector x) const
-        {
-            return static_cast<T>(rosen(x));
-        }
-
-        void get_derivative_and_hessian (
-            const column_vector& x,
-            column_vector& d,
-            general_matrix& h
-        ) const 
-        {
-            d = rosen_derivative(x);
-            h = rosen_hessian(x);
-        }
-
-    };
-
 // ----------------------------------------------------------------------------------------
 
     template <typename T>
@@ -131,7 +63,7 @@ namespace
 
         matrix<T,2,1> p = 100*matrix_cast<T>(randm(2,1,rnd)) - 50;
 
-        T obj = find_min_trust_region(objective_delta_stop_strategy(0, 100), rosen_model<T>(), p);
+        T obj = find_min_trust_region(objective_delta_stop_strategy(0, 100), rosen_function_model<T>(), p);
 
         DLIB_TEST_MSG(obj == 0, "obj: " << obj);
         DLIB_TEST_MSG(length(p-ans) == 0, "length(p): " << length(p-ans));
@@ -258,6 +190,113 @@ namespace
 
 // ----------------------------------------------------------------------------------------
 
+    void test_problems()
+    {
+        print_spinner();
+        {
+            matrix<double,4,1> ch;
+
+            ch = brown_start();
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-7, 80),
+                                  brown_function_model(),
+                                  ch);
+
+            dlog << LINFO << "brown obj: " << brown(ch);
+            dlog << LINFO << "brown der: " << length(brown_derivative(ch));
+            dlog << LINFO << "brown error: " << length(ch - brown_solution());
+
+            DLIB_TEST(length(ch - brown_solution()) < 1e-5);
+
+        }
+        print_spinner();
+        {
+            matrix<double,2,1> ch;
+
+            ch = rosen_start<double>();
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-7, 80),
+                                  rosen_function_model<double>(),
+                                  ch);
+
+            dlog << LINFO << "rosen obj: " << rosen(ch);
+            dlog << LINFO << "rosen der: " << length(rosen_derivative(ch));
+            dlog << LINFO << "rosen error: " << length(ch - rosen_solution<double>());
+
+            DLIB_TEST(length(ch - rosen_solution<double>()) < 1e-5);
+        }
+
+        print_spinner();
+        {
+            matrix<double,0,1> ch;
+
+            ch = chebyquad_start(2);
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-7, 80),
+                                  chebyquad_function_model(),
+                                  ch);
+
+            dlog << LINFO << "chebyquad 2 obj: " << chebyquad(ch);
+            dlog << LINFO << "chebyquad 2 der: " << length(chebyquad_derivative(ch));
+            dlog << LINFO << "chebyquad 2 error: " << length(ch - chebyquad_solution(2));
+
+            DLIB_TEST(length(ch - chebyquad_solution(2)) < 1e-5);
+
+        }
+        print_spinner();
+        {
+            matrix<double,0,1> ch;
+
+            ch = chebyquad_start(4);
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-7, 80),
+                                  chebyquad_function_model(),
+                                  ch);
+
+            dlog << LINFO << "chebyquad 4 obj: " << chebyquad(ch);
+            dlog << LINFO << "chebyquad 4 der: " << length(chebyquad_derivative(ch));
+            dlog << LINFO << "chebyquad 4 error: " << length(ch - chebyquad_solution(4));
+
+            DLIB_TEST(length(ch - chebyquad_solution(4)) < 1e-5);
+        }
+        print_spinner();
+        {
+            matrix<double,0,1> ch;
+
+            ch = chebyquad_start(6);
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-12, 80),
+                                  chebyquad_function_model(),
+                                  ch);
+
+            dlog << LINFO << "chebyquad 6 obj: " << chebyquad(ch);
+            dlog << LINFO << "chebyquad 6 der: " << length(chebyquad_derivative(ch));
+            dlog << LINFO << "chebyquad 6 error: " << length(ch - chebyquad_solution(6));
+
+            DLIB_TEST(length(ch - chebyquad_solution(6)) < 1e-5);
+
+        }
+        print_spinner();
+        {
+            matrix<double,0,1> ch;
+
+            ch = chebyquad_start(8);
+
+            find_min_trust_region(objective_delta_stop_strategy(1e-10, 80),
+                                  chebyquad_function_model(),
+                                  ch);
+
+            dlog << LINFO << "chebyquad 8 obj: " << chebyquad(ch);
+            dlog << LINFO << "chebyquad 8 der: " << length(chebyquad_derivative(ch));
+            dlog << LINFO << "chebyquad 8 error: " << length(ch - chebyquad_solution(8));
+
+            DLIB_TEST(length(ch - chebyquad_solution(8)) < 1e-5);
+        }
+
+    }
+
+
+
     class optimization_tester : public tester
     {
     public:
@@ -280,6 +319,8 @@ namespace
 
 
             test_trust_region_sub_problem();
+
+            test_problems();
         }
     } a;