diff --git a/fdaPDE/fields/scalar_field.h b/fdaPDE/fields/scalar_field.h index 08507a23..ff32d061 100644 --- a/fdaPDE/fields/scalar_field.h +++ b/fdaPDE/fields/scalar_field.h @@ -47,6 +47,7 @@ class ScalarField : public ScalarExpr> { typedef typename static_dynamic_vector_selector::type VectorType; typedef typename static_dynamic_matrix_selector::type MatrixType; typedef ScalarExpr> Base; + static constexpr int DomainDimension = N; static_assert( std::is_invocable::value && std::is_same::type, double>::value); @@ -114,6 +115,7 @@ template struct ScalarField : // non-const member function pointers public ScalarField_MemFnBase { using Base = ScalarField_MemFnBase; + static constexpr int DomainDimension = N; // constructors ScalarField() = default; ScalarField(typename Base::ClassPtrType_ c, typename Base::FieldType_ f) : Base(c, f) {}; diff --git a/fdaPDE/optimization.h b/fdaPDE/optimization.h index ae53eb12..e9e23ace 100644 --- a/fdaPDE/optimization.h +++ b/fdaPDE/optimization.h @@ -17,6 +17,7 @@ #ifndef __FDAPDE_OPTIMIZATION_MODULE_H__ #define __FDAPDE_OPTIMIZATION_MODULE_H__ +#include "optimization/optimizer.h" #include "optimization/grid.h" #include "optimization/newton.h" #include "optimization/gradient_descent.h" diff --git a/fdaPDE/optimization/bfgs.h b/fdaPDE/optimization/bfgs.h index 540ed34e..f3fdf2c8 100644 --- a/fdaPDE/optimization/bfgs.h +++ b/fdaPDE/optimization/bfgs.h @@ -25,10 +25,11 @@ namespace fdapde { namespace core { // implementation of the Broyden–Fletcher–Goldfarb–Shanno algorithm for unconstrained nonlinear optimization -template class BFGS { +template class BFGS { private: typedef typename std::conditional, SVector>::type VectorType; typedef typename std::conditional, SMatrix>::type MatrixType; + std::tuple callbacks_ {}; std::size_t max_iter_; // maximum number of iterations before forced stop double tol_; // tolerance on error before forced stop double step_; // update step @@ -43,11 +44,13 @@ template class BFGS { // constructor BFGS() = default; BFGS(std::size_t max_iter, double tol, double step) : max_iter_(max_iter), tol_(tol), step_(step) {}; + BFGS(std::size_t max_iter, double tol, double step, Args&... callbacks) : + max_iter_(max_iter), tol_(tol), step_(step), callbacks_(std::make_tuple(std::forward(callbacks)...)) {}; - template VectorType optimize(F& objective, const VectorType& x0, Args... args) { + template VectorType optimize(F& objective, const VectorType& x0) { static_assert( std::is_same().operator()(VectorType())), double>::value, - "cannot find definition for F.operator()(const VectorType&)"); + "F_IS_NOT_A_FUNCTOR_ACCEPTING_A_VECTORTYPE"); bool stop = false; // asserted true in case of forced stop VectorType zero; // either statically or dynamically allocated depending on N @@ -72,7 +75,7 @@ template class BFGS { while (n_iter < max_iter_ && error > tol_ && !stop) { // compute update direction update = -inv_hessian * grad_old; - stop |= execute_pre_update_step(*this, objective, args...); + stop |= execute_pre_update_step(*this, objective, callbacks_); // update along descent direction x_new = x_old + h * update; @@ -95,7 +98,7 @@ template class BFGS { // prepare next iteration error = grad_new.norm(); - stop |= execute_post_update_step(*this, objective, args...); + stop |= execute_post_update_step(*this, objective, callbacks_); x_old = x_new; grad_old = grad_new; n_iter++; diff --git a/fdaPDE/optimization/callbacks/callbacks.h b/fdaPDE/optimization/callbacks/callbacks.h index 9f99a4af..048586b7 100644 --- a/fdaPDE/optimization/callbacks/callbacks.h +++ b/fdaPDE/optimization/callbacks/callbacks.h @@ -27,28 +27,26 @@ define_has(pre_update_step); define_has(post_update_step); template -bool execute_pre_update_step(Opt& optimizer, Obj& objective, Args&... args) { +bool execute_pre_update_step(Opt& optimizer, Obj& objective, std::tuple& callbacks) { bool b = false; - ( // fold expand parameter pack - [&] { - if constexpr (has_pre_update_step::value) { - b |= args.pre_update_step(optimizer, objective); - } - }(), - ...); + auto exec_callback = [&](auto&& callback) { + if constexpr (has_pre_update_step, bool(Opt&, Obj&)>::value) { + b |= callback.pre_update_step(optimizer, objective); + } + }; + std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks); return b; } template -bool execute_post_update_step(Opt& optimizer, Obj& objective, Args&... args) { +bool execute_post_update_step(Opt& optimizer, Obj& objective, std::tuple& callbacks) { bool b = false; - ( // fold expand parameter pack - [&] { - if constexpr (has_post_update_step::value) { - b |= args.post_update_step(optimizer, objective); - } - }(), - ...); + auto exec_callback = [&](auto&& callback) { + if constexpr (has_post_update_step, bool(Opt&, Obj&)>::value) { + b |= callback.post_update_step(optimizer, objective); + } + }; + std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks); return b; } diff --git a/fdaPDE/optimization/callbacks/grid_start.h b/fdaPDE/optimization/callbacks/grid_start.h new file mode 100644 index 00000000..9efea08f --- /dev/null +++ b/fdaPDE/optimization/callbacks/grid_start.h @@ -0,0 +1,53 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __GRID_START_H__ +#define __GRID_START_H__ + +#include "../../utils/symbols.h" + +namespace fdapde { +namespace core { + +// implementation of the backatracking line search method for step selection +class GridStart { + private: + double alpha_ = 2.0; + double beta_ = 0.5; + double gamma_ = 0.5; + public: + // constructors + GridStart(double alpha, double beta, double gamma) : alpha_(alpha), beta_(beta), gamma_(gamma) {}; + + // backtracking based step search + template bool pre_update_step(Opt& opt, Obj& obj) { + double alpha = alpha_; // restore to user defined settings + double m = opt.grad_old.dot(opt.update); + if (m < 0) { // descent direction + while (obj(opt.x_old) - obj(opt.x_old + alpha * opt.update) // Armijo–Goldstein condition + + gamma_ * alpha * m < 0) { + alpha *= beta_; + } + } + opt.h = alpha; + return false; + } +}; + +} // namespace core +} // namespace fdapde + +#endif // __GRID_START_H__ diff --git a/fdaPDE/optimization/gradient_descent.h b/fdaPDE/optimization/gradient_descent.h index 6e53a14b..3dbbd3c7 100644 --- a/fdaPDE/optimization/gradient_descent.h +++ b/fdaPDE/optimization/gradient_descent.h @@ -25,10 +25,11 @@ namespace fdapde { namespace core { // implementation of the gradient descent method for unconstrained nonlinear optimization -template class GradientDescent { +template class GradientDescent { private: typedef typename std::conditional, SVector>::type VectorType; typedef typename std::conditional, SMatrix>::type MatrixType; + std::tuple callbacks_ {}; std::size_t max_iter_; // maximum number of iterations before forced stop double tol_; // tolerance on error before forced stop double step_; // update step @@ -43,11 +44,13 @@ template class GradientDescent { // constructor GradientDescent() = default; GradientDescent(std::size_t max_iter, double tol, double step) : max_iter_(max_iter), tol_(tol), step_(step) {}; + GradientDescent(std::size_t max_iter, double tol, double step, Args&... callbacks) : + max_iter_(max_iter), tol_(tol), step_(step), callbacks_(std::make_tuple(std::forward(callbacks)...)) {}; - template VectorType optimize(F& objective, const VectorType& x0, Args... args) { + template VectorType optimize(F& objective, const VectorType& x0) { static_assert( std::is_same().operator()(VectorType())), double>::value, - "cannot find definition for F.operator()(const VectorType&)"); + "F_IS_NOT_A_FUNCTOR_ACCEPTING_A_VECTORTYPE"); bool stop = false; // asserted true in case of forced stop std::size_t n_iter = 0; @@ -60,7 +63,7 @@ template class GradientDescent { while (n_iter < max_iter_ && error > tol_ && !stop) { update = -grad_old; - stop |= execute_pre_update_step(*this, objective, args...); + stop |= execute_pre_update_step(*this, objective, callbacks_); // update along descent direction x_new = x_old + h * update; @@ -68,7 +71,7 @@ template class GradientDescent { // prepare next iteration error = grad_new.norm(); - stop |= execute_post_update_step(*this, objective, args...); + stop |= execute_post_update_step(*this, objective, callbacks_); x_old = x_new; grad_old = grad_new; n_iter++; diff --git a/fdaPDE/optimization/grid.h b/fdaPDE/optimization/grid.h index b8eb3ab0..a31abd82 100644 --- a/fdaPDE/optimization/grid.h +++ b/fdaPDE/optimization/grid.h @@ -25,22 +25,24 @@ namespace fdapde { namespace core { // searches for the point in a given grid minimizing a given nonlinear objective -template class Grid { +template class Grid { private: typedef typename std::conditional, SVector>::type VectorType; + std::tuple callbacks_ {}; VectorType optimum_; double value_; // objective value at optimum public: VectorType x_current; // constructor - Grid() = default; + template ::type = 0> Grid() {}; + Grid(Args&&... callbacks) : callbacks_(std::make_tuple(std::forward(callbacks)...)) {}; - template - VectorType optimize(F& objective, const std::vector& grid, Args&... args) { + template + VectorType optimize(F& objective, const std::vector& grid) { static_assert( std::is_same().operator()(VectorType())), double>::value, - "cannot find definition for F.operator()(const VectorType&)"); + "F_IS_NOT_A_FUNCTOR_ACCEPTING_A_VECTORTYPE"); bool stop = false; // asserted true in case of forced stop // algorithm initialization @@ -51,7 +53,7 @@ template class Grid { for (std::size_t i = 1; i < grid.size() && !stop; ++i) { x_current = grid[i]; double x = objective(x_current); - stop |= execute_post_update_step(*this, objective, args...); + stop |= execute_post_update_step(*this, objective, callbacks_); // update minimum if better optimum found if (x < value_) { diff --git a/fdaPDE/optimization/newton.h b/fdaPDE/optimization/newton.h index a2b1a826..e5933e10 100644 --- a/fdaPDE/optimization/newton.h +++ b/fdaPDE/optimization/newton.h @@ -25,10 +25,11 @@ namespace fdapde { namespace core { // implementation of the newton method for unconstrained nonlinear optimization -template class Newton { +template class Newton { private: typedef typename std::conditional, SVector>::type VectorType; typedef typename std::conditional, SMatrix>::type MatrixType; + std::tuple callbacks_ {}; std::size_t max_iter_; // maximum number of iterations before forced stop double tol_; // tolerance on error before forced stop double step_; // update step @@ -44,11 +45,13 @@ template class Newton { // constructor Newton() = default; Newton(std::size_t max_iter, double tol, double step) : max_iter_(max_iter), tol_(tol), step_(step) {}; + Newton(std::size_t max_iter, double tol, double step, Args&... callbacks) : + max_iter_(max_iter), tol_(tol), step_(step), callbacks_(std::make_tuple(std::forward(callbacks)...)) {}; - template VectorType optimize(F& objective, const VectorType& x0, Args... args) { + template VectorType optimize(F& objective, const VectorType& x0) { static_assert( std::is_same().operator()(VectorType())), double>::value, - "cannot find definition for F.operator()(const VectorType&)"); + "F_IS_NOT_A_FUNCTOR_ACCEPTING_A_VECTORTYPE"); bool stop = false; // asserted true in case of forced stop std::size_t n_iter = 0; @@ -64,7 +67,7 @@ template class Newton { inv_hessian.compute(hessian); update = -inv_hessian.solve(grad_old); - stop |= execute_pre_update_step(*this, objective, args...); + stop |= execute_pre_update_step(*this, objective, callbacks_); // update along descent direction x_new = x_old + h * update; @@ -72,7 +75,7 @@ template class Newton { // prepare next iteration error = grad_new.norm(); - stop |= execute_post_update_step(*this, objective, args...); + stop |= execute_post_update_step(*this, objective, callbacks_); x_old = x_new; grad_old = grad_new; n_iter++; diff --git a/fdaPDE/optimization/optimizer.h b/fdaPDE/optimization/optimizer.h new file mode 100644 index 00000000..efcaea98 --- /dev/null +++ b/fdaPDE/optimization/optimizer.h @@ -0,0 +1,44 @@ +// This file is part of fdaPDE, a C++ library for physics-informed +// spatial and functional data analysis. +// +// This program is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program 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 General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef __OPTIMIZER_H__ +#define __OPTIMIZER_H__ + +#include "../utils/symbols.h" +#include "../utils/type_erasure.h" + +namespace fdapde { +namespace core { + +// a type-erasure wrapper for an optimization algorithm optimizing objectives of type F +template struct Optimizer__ { + static constexpr int N = F::DomainDimension; + using VectorType = typename std::conditional, SVector>::type; + // function pointers forwardings + template using fn_ptrs = fdapde::mem_fn_ptrs<&O::template optimize, &O::optimum, &O::value>; + // implementation + template VectorType optimize(F& objective, const T& x0) { + return fdapde::invoke(*this, objective, x0); + } + VectorType optimum() { return fdapde::invoke(*this); } + double value() { return fdapde::invoke(*this); } +}; +template using Optimizer = fdapde::erase>; + +} // namespace core +} // namespace fdapde + +#endif // __OPTIMIZER_H__ diff --git a/test/src/optimization_test.cpp b/test/src/optimization_test.cpp index 68fd581a..7a54ff94 100644 --- a/test/src/optimization_test.cpp +++ b/test/src/optimization_test.cpp @@ -20,6 +20,7 @@ #include #include #include +using fdapde::core::Optimizer; using fdapde::core::BFGS; using fdapde::core::GradientDescent; using fdapde::core::Grid; @@ -57,9 +58,9 @@ TEST(optimization_test, gradient_descent_backtracking_line_search) { }; f.set_step(1e-4); - GradientDescent<2> opt(1000, 1e-6, 0.01); + GradientDescent<2, BacktrackingLineSearch> opt(1000, 1e-6, 0.01); SVector<2> pt(-1, -1); - opt.optimize(f, pt, BacktrackingLineSearch()); + opt.optimize(f, pt); // expected solution SVector<2> expected(-0.6690718221499544, 0); @@ -75,9 +76,9 @@ TEST(optimization_test, newton_backtracking_line_search) { }; f.set_step(1e-4); - Newton<2> opt(1000, 1e-6, 0.01); + Newton<2, BacktrackingLineSearch> opt(1000, 1e-6, 0.01); SVector<2> pt(-0.5, -0.5); - opt.optimize(f, pt, BacktrackingLineSearch()); + opt.optimize(f, pt); // expected solution SVector<2> expected(-0.6690718221499544, 0); @@ -86,7 +87,7 @@ TEST(optimization_test, newton_backtracking_line_search) { } // test integration of Laplacian weak form for a LagrangianBasis of order 2 -TEST(optimization_test, bfgs_wolfe_line_search) { +TEST(optimization_test, type_erased_bfgs_wolfe_line_search) { // define objective function: x*e^{-x^2 - y^2} + (x^2 + y^2)/20 ScalarField<2> f; f = [](SVector<2> x) -> double { @@ -95,10 +96,10 @@ TEST(optimization_test, bfgs_wolfe_line_search) { f.set_step(1e-4); // define optimizer - BFGS<2> opt(1000, 1e-6, 0.01); + Optimizer> opt = BFGS<2, WolfeLineSearch>(1000, 1e-6, 0.01); // use a type erasure wrapper // perform optimization SVector<2> pt(-1, -1); - opt.optimize(f, pt, WolfeLineSearch()); + opt.optimize(f, pt); // expected solution SVector<2> expected(-0.6690718221499544, 0);