Skip to content

Commit

Permalink
optimizer's callbacks passed via constructor, (now are a parameter pa…
Browse files Browse the repository at this point in the history
…ck of an optimizer template, stored inside the class), type_erasure wrapper for a general concept of optimizer, scalar_field exposes its domain dimensionality at compile time
  • Loading branch information
AlePalu committed Jan 19, 2024
1 parent 182556d commit 93ce5f8
Show file tree
Hide file tree
Showing 10 changed files with 154 additions and 44 deletions.
2 changes: 2 additions & 0 deletions fdaPDE/fields/scalar_field.h
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ class ScalarField : public ScalarExpr<N, ScalarField<N, F>> {
typedef typename static_dynamic_vector_selector<N>::type VectorType;
typedef typename static_dynamic_matrix_selector<N, N>::type MatrixType;
typedef ScalarExpr<N, ScalarField<N, F>> Base;
static constexpr int DomainDimension = N;
static_assert(
std::is_invocable<F, VectorType>::value &&
std::is_same<typename std::invoke_result<F, VectorType>::type, double>::value);
Expand Down Expand Up @@ -114,6 +115,7 @@ template <int N, typename RetType_, typename ClassType_, typename VectorType_>
struct ScalarField<N, RetType_ (ClassType_::*)(VectorType_)> : // non-const member function pointers
public ScalarField_MemFnBase<N, RetType_ (ClassType_::*)(VectorType_)> {
using Base = ScalarField_MemFnBase<N, RetType_ (ClassType_::*)(VectorType_)>;
static constexpr int DomainDimension = N;
// constructors
ScalarField() = default;
ScalarField(typename Base::ClassPtrType_ c, typename Base::FieldType_ f) : Base(c, f) {};
Expand Down
1 change: 1 addition & 0 deletions fdaPDE/optimization.h
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
13 changes: 8 additions & 5 deletions fdaPDE/optimization/bfgs.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,11 @@ namespace fdapde {
namespace core {

// implementation of the Broyden–Fletcher–Goldfarb–Shanno algorithm for unconstrained nonlinear optimization
template <int N> class BFGS {
template <int N, typename... Args> class BFGS {
private:
typedef typename std::conditional<N == Dynamic, DVector<double>, SVector<N>>::type VectorType;
typedef typename std::conditional<N == Dynamic, DMatrix<double>, SMatrix<N>>::type MatrixType;
std::tuple<Args...> 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
Expand All @@ -43,11 +44,13 @@ template <int N> 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<Args>(callbacks)...)) {};

template <typename F, typename... Args> VectorType optimize(F& objective, const VectorType& x0, Args... args) {
template <typename F> VectorType optimize(F& objective, const VectorType& x0) {
static_assert(
std::is_same<decltype(std::declval<F>().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
Expand All @@ -72,7 +75,7 @@ template <int N> 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;
Expand All @@ -95,7 +98,7 @@ template <int N> 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++;
Expand Down
30 changes: 14 additions & 16 deletions fdaPDE/optimization/callbacks/callbacks.h
Original file line number Diff line number Diff line change
Expand Up @@ -27,28 +27,26 @@ define_has(pre_update_step);
define_has(post_update_step);

template <typename Opt, typename Obj, typename... Args>
bool execute_pre_update_step(Opt& optimizer, Obj& objective, Args&... args) {
bool execute_pre_update_step(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
bool b = false;
( // fold expand parameter pack
[&] {
if constexpr (has_pre_update_step<Args, bool(Opt&, Obj&)>::value) {
b |= args.pre_update_step(optimizer, objective);
}
}(),
...);
auto exec_callback = [&](auto&& callback) {
if constexpr (has_pre_update_step<std::decay_t<decltype(callback)>, bool(Opt&, Obj&)>::value) {
b |= callback.pre_update_step(optimizer, objective);
}
};
std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks);
return b;
}

template <typename Opt, typename Obj, typename... Args>
bool execute_post_update_step(Opt& optimizer, Obj& objective, Args&... args) {
bool execute_post_update_step(Opt& optimizer, Obj& objective, std::tuple<Args...>& callbacks) {
bool b = false;
( // fold expand parameter pack
[&] {
if constexpr (has_post_update_step<Args, bool(Opt&, Obj&)>::value) {
b |= args.post_update_step(optimizer, objective);
}
}(),
...);
auto exec_callback = [&](auto&& callback) {
if constexpr (has_post_update_step<std::decay_t<decltype(callback)>, bool(Opt&, Obj&)>::value) {
b |= callback.post_update_step(optimizer, objective);
}
};
std::apply([&](auto&&... callback) { (exec_callback(callback), ...); }, callbacks);
return b;
}

Expand Down
53 changes: 53 additions & 0 deletions fdaPDE/optimization/callbacks/grid_start.h
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.

#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 <typename Opt, typename Obj> 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__
13 changes: 8 additions & 5 deletions fdaPDE/optimization/gradient_descent.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,11 @@ namespace fdapde {
namespace core {

// implementation of the gradient descent method for unconstrained nonlinear optimization
template <int N> class GradientDescent {
template <int N, typename... Args> class GradientDescent {
private:
typedef typename std::conditional<N == Dynamic, DVector<double>, SVector<N>>::type VectorType;
typedef typename std::conditional<N == Dynamic, DMatrix<double>, SMatrix<N>>::type MatrixType;
std::tuple<Args...> 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
Expand All @@ -43,11 +44,13 @@ template <int N> 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<Args>(callbacks)...)) {};

template <typename F, typename... Args> VectorType optimize(F& objective, const VectorType& x0, Args... args) {
template <typename F> VectorType optimize(F& objective, const VectorType& x0) {
static_assert(
std::is_same<decltype(std::declval<F>().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;
Expand All @@ -60,15 +63,15 @@ template <int N> 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;
grad_new = objective.derive()(x_new);

// 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++;
Expand Down
14 changes: 8 additions & 6 deletions fdaPDE/optimization/grid.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,22 +25,24 @@ namespace fdapde {
namespace core {

// searches for the point in a given grid minimizing a given nonlinear objective
template <int N> class Grid {
template <int N, typename... Args> class Grid {
private:
typedef typename std::conditional<N == Dynamic, DVector<double>, SVector<N>>::type VectorType;
std::tuple<Args...> callbacks_ {};
VectorType optimum_;
double value_; // objective value at optimum
public:
VectorType x_current;

// constructor
Grid() = default;
template <int N_ = sizeof...(Args), typename std::enable_if<N_ != 0, int>::type = 0> Grid() {};
Grid(Args&&... callbacks) : callbacks_(std::make_tuple(std::forward<Args>(callbacks)...)) {};

template <typename F, typename... Args>
VectorType optimize(F& objective, const std::vector<VectorType>& grid, Args&... args) {
template <typename F>
VectorType optimize(F& objective, const std::vector<VectorType>& grid) {
static_assert(
std::is_same<decltype(std::declval<F>().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
Expand All @@ -51,7 +53,7 @@ template <int N> 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_) {
Expand Down
13 changes: 8 additions & 5 deletions fdaPDE/optimization/newton.h
Original file line number Diff line number Diff line change
Expand Up @@ -25,10 +25,11 @@ namespace fdapde {
namespace core {

// implementation of the newton method for unconstrained nonlinear optimization
template <int N> class Newton {
template <int N, typename... Args> class Newton {
private:
typedef typename std::conditional<N == Dynamic, DVector<double>, SVector<N>>::type VectorType;
typedef typename std::conditional<N == Dynamic, DMatrix<double>, SMatrix<N>>::type MatrixType;
std::tuple<Args...> 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
Expand All @@ -44,11 +45,13 @@ template <int N> 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<Args>(callbacks)...)) {};

template <typename F, typename... Args> VectorType optimize(F& objective, const VectorType& x0, Args... args) {
template <typename F> VectorType optimize(F& objective, const VectorType& x0) {
static_assert(
std::is_same<decltype(std::declval<F>().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;
Expand All @@ -64,15 +67,15 @@ template <int N> 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;
grad_new = objective.derive()(x_new);

// 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++;
Expand Down
44 changes: 44 additions & 0 deletions fdaPDE/optimization/optimizer.h
Original file line number Diff line number Diff line change
@@ -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 <http://www.gnu.org/licenses/>.

#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 <typename F> struct Optimizer__ {
static constexpr int N = F::DomainDimension;
using VectorType = typename std::conditional<N == Dynamic, DVector<double>, SVector<N>>::type;
// function pointers forwardings
template <typename O> using fn_ptrs = fdapde::mem_fn_ptrs<&O::template optimize<F>, &O::optimum, &O::value>;
// implementation
template <typename T> VectorType optimize(F& objective, const T& x0) {
return fdapde::invoke<VectorType, 0>(*this, objective, x0);
}
VectorType optimum() { return fdapde::invoke<VectorType, 1>(*this); }
double value() { return fdapde::invoke<double, 2>(*this); }
};
template <typename F> using Optimizer = fdapde::erase<fdapde::heap_storage, Optimizer__<F>>;

} // namespace core
} // namespace fdapde

#endif // __OPTIMIZER_H__
15 changes: 8 additions & 7 deletions test/src/optimization_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <fdaPDE/utils.h>
#include <fdaPDE/fields.h>
#include <fdaPDE/optimization.h>
using fdapde::core::Optimizer;
using fdapde::core::BFGS;
using fdapde::core::GradientDescent;
using fdapde::core::Grid;
Expand Down Expand Up @@ -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);
Expand All @@ -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);
Expand All @@ -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 {
Expand All @@ -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<ScalarField<2>> 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);
Expand Down

0 comments on commit 93ce5f8

Please sign in to comment.