Example shows numeric quadratic function use.
// Copyright (C) 2009 by Thomas Moulard, AIST, CNRS, INRIA. // // This file is part of the roboptim. // // roboptim is free software: you can redistribute it and/or modify // it under the terms of the GNU Lesser General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. // // roboptim is distributed in the hope that it will be useful, // but WITHOUT ANY WARRANTY; without even the implied warranty of // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the // GNU Lesser General Public License for more details. // // You should have received a copy of the GNU Lesser General Public License // along with roboptim. If not, see <http://www.gnu.org/licenses/>. #include "shared-tests/fixture.hh" #include <iostream> #include <boost/type_traits/is_same.hpp> #include <roboptim/core/io.hh> #include <roboptim/core/decorator/finite-difference-gradient.hh> #include <roboptim/core/numeric-quadratic-function.hh> #include <roboptim/core/util.hh> using namespace roboptim; typedef DummySolver solver_t; typedef boost::mpl::list< ::roboptim::EigenMatrixDense, ::roboptim::EigenMatrixSparse> functionTypes_t; BOOST_FIXTURE_TEST_SUITE (core, TestSuiteConfiguration) BOOST_AUTO_TEST_CASE_TEMPLATE (numeric_quadratic_function, T, functionTypes_t) { boost::shared_ptr<boost::test_tools::output_test_stream> output = retrievePattern ("numeric-quadratic-function"); typename GenericNumericQuadraticFunction<T>::matrix_t a (5, 5); typename GenericNumericQuadraticFunction<T>::vector_t b (5); typename GenericNumericQuadraticFunction<T>::vector_t x (5); a.setZero (); b.setZero (); for (int i = 0; i < 5; ++i) { a.coeffRef (i,i) = 1.; x[i] = static_cast<Function::value_type> (i); } GenericNumericQuadraticFunction<T> f (a, b, "Dummy"); (*output) << f << '\n'; (*output) << "x = " << x << '\n'; (*output) << "f(x) = " << f (x) << '\n'; (*output) << "J(x) = " << toDense (f.jacobian (x)) << '\n'; (*output) << "G(x) = " << toDense (f.gradient (x, 0)) << '\n'; (*output) << "H(x) = " << toDense (f.hessian (x, 0)) << '\n'; std::cout << output->str () << std::endl; if (boost::is_same<T, EigenMatrixDense>::value) BOOST_CHECK (output->match_pattern ()); for (int i = 0; i < 10; ++i) { for (int j = 0; j < 5; ++j) x[j] = std::ceil (rand () % 50); std::cout << "x = " << x << '\n'; std::cout << "f(x) = " << f (x) << '\n'; std::cout << "J(x) = " << f.jacobian (x) << '\n'; std::cout << "G(x) = " << f.gradient (x, 0) << '\n'; std::cout << "H(x) = " << f.hessian (x, 0) << '\n'; typedef typename GenericNumericQuadraticFunction<T>::matrix_t matrix_t; matrix_t J (1, 5); for (typename matrix_t::Index i = 0; i < 5; ++i) J.coeffRef (0, i) = 2 * x[i]; for (typename matrix_t::Index i = 0; i < 5; ++i) { std::cout << f.jacobian (x).coeffRef (0, i) << '\n'; std::cout << J.coeffRef (0, i) << '\n'; } BOOST_CHECK (allclose (f.jacobian (x), J)); BOOST_CHECK (allclose (f.hessian (x, 0), 2*a)); BOOST_CHECK (checkGradient (f, 0, x)); BOOST_CHECK (checkJacobian (f, x)); } } typedef boost::mpl::list< ::roboptim::EigenMatrixSparse> sparseOnly_t; BOOST_AUTO_TEST_CASE_TEMPLATE (random_gradient_check, T, sparseOnly_t) { typename GenericNumericQuadraticFunction<T>::matrix_t a (5, 5); typename GenericNumericQuadraticFunction<T>::vector_t b (5); typename GenericNumericQuadraticFunction<T>::vector_t x (5); for (int randomTry = 0; randomTry < 10; ++randomTry) { a.setZero (); b.setZero (); x.setZero (); for (typename GenericNumericQuadraticFunction<T>::matrix_t::Index i = 0; i < 5; ++i) for (typename GenericNumericQuadraticFunction<T>::matrix_t::Index j = 0; j < 5; ++j) a.insert (i, j) = 0.; for (typename GenericNumericQuadraticFunction<T>::matrix_t::Index i = 0; i < 5; ++i) for (typename GenericNumericQuadraticFunction<T>::matrix_t::Index j = 0; j < 5; ++j) a.coeffRef (i, j) = a.coeffRef (j, i) = static_cast<double> (std::rand () / RAND_MAX); b = GenericNumericQuadraticFunction<T>::vector_t::Random (5); GenericNumericQuadraticFunction<T> f (a, b); for (int i = 0; i < 10; ++i) { for (int j = 0; j < 5; ++j) x[j] = std::ceil (rand () % 50); BOOST_CHECK (checkGradient (f, 0, x)); BOOST_CHECK (checkJacobian (f, x)); } } } BOOST_AUTO_TEST_SUITE_END ()