From 64470d480bd6b89c08356af832a1e0177b6ad70e Mon Sep 17 00:00:00 2001 From: imraneamri Date: Wed, 5 Jun 2024 15:46:27 +0100 Subject: [PATCH 01/30] Groundwork for Hessian class --- src/XAD/Hessian.hpp | 80 +++++++++++++++++++++++++++++++++++++++++++ src/XAD/XAD.hpp | 1 + test/CMakeLists.txt | 1 + test/Hessian_test.cpp | 59 +++++++++++++++++++++++++++++++ 4 files changed, 141 insertions(+) create mode 100644 src/XAD/Hessian.hpp create mode 100644 test/Hessian_test.cpp diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp new file mode 100644 index 0000000..e403622 --- /dev/null +++ b/src/XAD/Hessian.hpp @@ -0,0 +1,80 @@ +/******************************************************************************* + + Implementation of Hessian computing methods. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +namespace xad +{ + +template +class Hessian +{ + public: + Hessian() {} + + std::vector> compute(std::vector &v) + { + xad::Tape> tape; + + domain = static_cast(v.size()); + std::vector> matrix(domain, std::vector(domain, 0.0)); + + tape.registerInputs(v); + + for (unsigned int i = 0; i < domain; i++) + { + tape.registerInputs(v); + derivative(value(v[i])) = 1.0; + tape.newRecording(); + + T y = q(v[0], v[1]); + tape.registerOutput(y); + value(derivative(y)) = 1.0; + + tape.computeAdjoints(); + + for (unsigned int j = 0; j < domain; j++) + { + std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) + << "\n"; + matrix[i][j] = derivative(derivative(v[j])); + } + + derivative(value(v[i])) = 0.0; + } + + return matrix; + } + + private: + T q(T a, T b) + { + T c = a * a; + T d = b * b; + return c + d; + } + + unsigned int domain; +}; +} // namespace xad diff --git a/src/XAD/XAD.hpp b/src/XAD/XAD.hpp index 0a1a310..088d1ed 100644 --- a/src/XAD/XAD.hpp +++ b/src/XAD/XAD.hpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include #include diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index de432ac..230e8cf 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -44,6 +44,7 @@ set(testfiles StdCompatibility_test.cpp ReusableRange_test.cpp PartialRollback_test.cpp + Hessian_test.cpp ) if (CMAKE_COMPILER_IS_GNUCC) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp new file mode 100644 index 0000000..2dd6635 --- /dev/null +++ b/test/Hessian_test.cpp @@ -0,0 +1,59 @@ +/******************************************************************************* + + Tests for xad::Hessian class. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +#include +#include +#include +#include +#include +#include + +using namespace ::testing; + +template +T quad(T a, T b) +{ + T c = a * a; + T d = b * b; + return c + d; +} + +TEST(HessianTest, SimpleQuadratic) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + + std::vector x = {3, 2}; + + xad::Hessian h; + + std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::vector> computed_hessian = h.compute(x); + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} From b0f8a7739821358051942c242536ac52e5a5e161 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Wed, 5 Jun 2024 17:20:21 +0100 Subject: [PATCH 02/30] Flexible function handling for hessian class --- src/XAD/Hessian.hpp | 16 +++++----------- test/Hessian_test.cpp | 35 ++++++++++++++++++++++++++++++----- 2 files changed, 35 insertions(+), 16 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index e403622..7721feb 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -31,7 +31,7 @@ template class Hessian { public: - Hessian() {} + Hessian(std::function &)> func) : foo(func) {} std::vector> compute(std::vector &v) { @@ -48,7 +48,7 @@ class Hessian derivative(value(v[i])) = 1.0; tape.newRecording(); - T y = q(v[0], v[1]); + T y = foo(v); tape.registerOutput(y); value(derivative(y)) = 1.0; @@ -56,8 +56,8 @@ class Hessian for (unsigned int j = 0; j < domain; j++) { - std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) - << "\n"; + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) + // << "\n"; matrix[i][j] = derivative(derivative(v[j])); } @@ -68,13 +68,7 @@ class Hessian } private: - T q(T a, T b) - { - T c = a * a; - T d = b * b; - return c + d; - } - + std::function &)> foo; unsigned int domain; }; } // namespace xad diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 2dd6635..2c103bd 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -34,21 +34,29 @@ using namespace ::testing; template -T quad(T a, T b) +T quad(std::vector &x) { - T c = a * a; - T d = b * b; + T c = x[0] * x[0]; + T d = x[1] * x[1]; return c + d; } +template +T tquad(std::vector &x) +{ + T c = x[0] * x[0]; + T d = x[1] * x[1]; + T e = x[2] * x[2]; + return c + d + e; +} + TEST(HessianTest, SimpleQuadratic) { typedef xad::fwd_adj mode; typedef mode::active_type AD; std::vector x = {3, 2}; - - xad::Hessian h; + xad::Hessian h(quad); std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; std::vector> computed_hessian = h.compute(x); @@ -57,3 +65,20 @@ TEST(HessianTest, SimpleQuadratic) for (unsigned int j = 0; j < cross_hessian.size(); j++) ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); } + +TEST(HessianTest, DDD) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + + std::vector x = {3, 2, 4}; + xad::Hessian h(tquad); + + std::vector> cross_hessian = { + {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}}; + std::vector> computed_hessian = h.compute(x); + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} From 22f3cc3eccc867ce86975bb1881aea1c499232e1 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Thu, 6 Jun 2024 10:30:09 +0100 Subject: [PATCH 03/30] Cleaned tests --- test/Hessian_test.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 2c103bd..77e6415 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -50,7 +50,16 @@ T tquad(std::vector &x) return c + d + e; } -TEST(HessianTest, SimpleQuadratic) +template +T foo(std::vector &x) +{ + T c = exp(x[0]); + T d = sin(x[1]); + T e = cos(x[2]); + return c + d + e; +} + +TEST(HessianTest, QuadraticForwardAdjoint) { typedef xad::fwd_adj mode; typedef mode::active_type AD; @@ -66,7 +75,7 @@ TEST(HessianTest, SimpleQuadratic) ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); } -TEST(HessianTest, DDD) +TEST(HessianTest, ThreeVarQuadratic) { typedef xad::fwd_adj mode; typedef mode::active_type AD; From 5807a1575af260f6dee7172f2413d1831eb54e10 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Thu, 6 Jun 2024 10:31:53 +0100 Subject: [PATCH 04/30] Added Hessian.hpp to CMakeList --- src/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 8f4b5ca..84fdc40 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -41,6 +41,7 @@ set(public_headers XAD/Complex.hpp XAD/Exceptions.hpp XAD/Expression.hpp + XAD/Hessian.hpp XAD/Interface.hpp XAD/Literals.hpp XAD/Macros.hpp From 33fef4f401a353d5c363186335778eb0f67acca2 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Thu, 6 Jun 2024 13:04:53 +0100 Subject: [PATCH 05/30] Jacobian base implementation --- src/CMakeLists.txt | 1 + src/XAD/Hessian.hpp | 8 ++--- src/XAD/Jacobian.hpp | 79 ++++++++++++++++++++++++++++++++++++++++++ src/XAD/XAD.hpp | 1 + test/CMakeLists.txt | 1 + test/Hessian_test.cpp | 11 +++--- test/Jacobian_test.cpp | 65 ++++++++++++++++++++++++++++++++++ 7 files changed, 156 insertions(+), 10 deletions(-) create mode 100644 src/XAD/Jacobian.hpp create mode 100644 test/Jacobian_test.cpp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 84fdc40..7b5321e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -43,6 +43,7 @@ set(public_headers XAD/Expression.hpp XAD/Hessian.hpp XAD/Interface.hpp + XAD/Jacobian.hpp XAD/Literals.hpp XAD/Macros.hpp XAD/MathFunctions.hpp diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 7721feb..878d628 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -22,6 +22,7 @@ ******************************************************************************/ +#include #include namespace xad @@ -31,12 +32,11 @@ template class Hessian { public: - Hessian(std::function &)> func) : foo(func) {} + Hessian(std::function &)> func, std::vector &v) : foo(func), v(v) {} - std::vector> compute(std::vector &v) + std::vector> compute() { xad::Tape> tape; - domain = static_cast(v.size()); std::vector> matrix(domain, std::vector(domain, 0.0)); @@ -44,7 +44,6 @@ class Hessian for (unsigned int i = 0; i < domain; i++) { - tape.registerInputs(v); derivative(value(v[i])) = 1.0; tape.newRecording(); @@ -69,6 +68,7 @@ class Hessian private: std::function &)> foo; + std::vector v; unsigned int domain; }; } // namespace xad diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp new file mode 100644 index 0000000..f0437ef --- /dev/null +++ b/src/XAD/Jacobian.hpp @@ -0,0 +1,79 @@ +/******************************************************************************* + + Implementation of Jacobian computing methods. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include +#include + +namespace xad +{ + +template +class Jacobian +{ + public: + Jacobian(std::vector &)>> foos, std::vector &v, + xad::Tape *tape) + : foos(foos), + v(v), + tape(tape), + domain(static_cast(v.size())), + codomain(static_cast(foos.size())) + { + } + + std::vector> compute() + { + std::vector> matrix(codomain, std::vector(domain, 0.0)); + + tape->registerInputs(v); + + for (unsigned int i = 0; i < domain; i++) + { + for (unsigned int j = 0; j < codomain; j++) + { + derivative(v[i]) = 1.0; + tape->newRecording(); + + T y = foos[j](v); + tape->registerOutput(y); + derivative(y) = 1.0; + + tape->computeAdjoints(); + + // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; + matrix[i][j] = derivative(v[i]); + derivative(v[i]) = 0.0; + } + } + + return matrix; + } + + private: + std::vector &)>> foos; + std::vector v; + Tape *tape; + unsigned int domain, codomain; +}; +} // namespace xad diff --git a/src/XAD/XAD.hpp b/src/XAD/XAD.hpp index 088d1ed..db912f0 100644 --- a/src/XAD/XAD.hpp +++ b/src/XAD/XAD.hpp @@ -37,6 +37,7 @@ #include #include #include +#include #include #include #include diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 230e8cf..0a03df3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,6 +45,7 @@ set(testfiles ReusableRange_test.cpp PartialRollback_test.cpp Hessian_test.cpp + Jacobian_test.cpp ) if (CMAKE_COMPILER_IS_GNUCC) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 77e6415..3dc54ed 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -61,14 +61,13 @@ T foo(std::vector &x) TEST(HessianTest, QuadraticForwardAdjoint) { - typedef xad::fwd_adj mode; - typedef mode::active_type AD; + typedef xad::AReal> AD; std::vector x = {3, 2}; - xad::Hessian h(quad); + xad::Hessian h(quad, x); std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> computed_hessian = h.compute(x); + std::vector> computed_hessian = h.compute(); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) @@ -81,11 +80,11 @@ TEST(HessianTest, ThreeVarQuadratic) typedef mode::active_type AD; std::vector x = {3, 2, 4}; - xad::Hessian h(tquad); + xad::Hessian hes(tquad, x); std::vector> cross_hessian = { {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}}; - std::vector> computed_hessian = h.compute(x); + std::vector> computed_hessian = hes.compute(); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp new file mode 100644 index 0000000..ad75d06 --- /dev/null +++ b/test/Jacobian_test.cpp @@ -0,0 +1,65 @@ +/******************************************************************************* + + Tests for xad::Jacobian class. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +#include +#include +#include +#include + +template +T foo1(std::vector &x) +{ + return x[0] + sin(x[1]); +} + +template +T foo2(std::vector &x) +{ + return x[1] + sin(x[0]); +} + +TEST(JacobianTest, SimpleJacobianAdjoint) +{ + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; + + std::vector &)>> funcs(2); + funcs[0] = foo1; + funcs[1] = foo2; + + std::vector x = {-2, 1}; + xad::Jacobian jac(funcs, x, &tape); + + std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; + std::vector> computed_jacobian = jac.compute(); + + for (unsigned int i = 0; i < cross_jacobian.size(); i++) + for (unsigned int j = 0; j < cross_jacobian.size(); j++) + ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); +} From 8dee7e85a61d22b6ff4b6afd068091ba32760756 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Thu, 6 Jun 2024 15:56:13 +0100 Subject: [PATCH 06/30] Added fwd_fwd mode for Hessian, changed Jacobian args + some new tests --- src/XAD/Hessian.hpp | 79 ++++++++++++++++++++++++++++++++---------- src/XAD/Jacobian.hpp | 56 +++++++++++++++++------------- test/Hessian_test.cpp | 62 ++++++++++++++++++++++++++++----- test/Jacobian_test.cpp | 14 ++++---- 4 files changed, 152 insertions(+), 59 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 878d628..6918941 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -32,43 +32,84 @@ template class Hessian { public: - Hessian(std::function &)> func, std::vector &v) : foo(func), v(v) {} + // fwd_adj constructor + Hessian(std::function &)> func, const std::vector &v, + xad::Tape> *tape) + : foo_(func), + v_(v), + domain_(static_cast(v_.size())), + matrix_(domain_, std::vector(domain_, 0.0)) + { + compute(tape); + } - std::vector> compute() + // tapeless (fwd_fwd) constructor + Hessian(std::function &)> func, std::vector &v) + : foo_(func), + v_(v), + domain_(static_cast(v_.size())), + matrix_(domain_, std::vector(domain_, 0.0)) { - xad::Tape> tape; - domain = static_cast(v.size()); - std::vector> matrix(domain, std::vector(domain, 0.0)); + compute(); + } - tape.registerInputs(v); + // fwd_adj + void compute(xad::Tape> *tape) + { + tape->registerInputs(v_); - for (unsigned int i = 0; i < domain; i++) + for (unsigned int i = 0; i < domain_; i++) { - derivative(value(v[i])) = 1.0; - tape.newRecording(); + derivative(value(v_[i])) = 1.0; + tape->newRecording(); - T y = foo(v); - tape.registerOutput(y); + T y = foo_(v_); + tape->registerOutput(y); value(derivative(y)) = 1.0; - tape.computeAdjoints(); + tape->computeAdjoints(); - for (unsigned int j = 0; j < domain; j++) + for (unsigned int j = 0; j < domain_; j++) { // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) // << "\n"; - matrix[i][j] = derivative(derivative(v[j])); + matrix_[i][j] = derivative(derivative(v_[j])); } - derivative(value(v[i])) = 0.0; + derivative(value(v_[i])) = 0.0; } + } - return matrix; + // fwd_fwd + void compute() + { + for (unsigned int i = 0; i < domain_; i++) + { + value(derivative(v_[i])) = 1.0; + + for (unsigned int j = 0; j < domain_; j++) + { + derivative(value(v_[j])) = 1.0; + + T y = foo_(v_); + + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) + // << "\n"; + matrix_[i][j] = derivative(derivative(y)); + + derivative(value(v_[j])) = 0.0; + } + + value(derivative(v_[i])) = 0.0; + } } + std::vector> get() { return matrix_; } + private: - std::function &)> foo; - std::vector v; - unsigned int domain; + std::function &)> foo_; + std::vector v_; + unsigned int domain_; + std::vector> matrix_; }; } // namespace xad diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index f0437ef..fffed58 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -32,48 +32,54 @@ template class Jacobian { public: - Jacobian(std::vector &)>> foos, std::vector &v, + Jacobian(std::function(std::vector &)> foo, const std::vector &v, xad::Tape *tape) - : foos(foos), - v(v), - tape(tape), - domain(static_cast(v.size())), - codomain(static_cast(foos.size())) + : foo_(foo), + v_(v), + tape_(tape), + domain_(static_cast(v_.size())), + codomain_(0), + matrix_(0) { + compute(); } - std::vector> compute() + void compute() { - std::vector> matrix(codomain, std::vector(domain, 0.0)); + tape_->registerInputs(v_); + std::vector res = foo_(v_); + codomain_ = static_cast(res.size()); + matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); - tape->registerInputs(v); - - for (unsigned int i = 0; i < domain; i++) + for (unsigned int i = 0; i < domain_; i++) { - for (unsigned int j = 0; j < codomain; j++) + for (unsigned int j = 0; j < codomain_; j++) { - derivative(v[i]) = 1.0; - tape->newRecording(); + derivative(v_[i]) = 1.0; + tape_->newRecording(); - T y = foos[j](v); - tape->registerOutput(y); + T y = res[j]; + tape_->registerOutput(y); derivative(y) = 1.0; - tape->computeAdjoints(); + tape_->computeAdjoints(); - // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; - matrix[i][j] = derivative(v[i]); - derivative(v[i]) = 0.0; + std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; + matrix_[i][j] = derivative(v_[i]); + derivative(v_[i]) = 0.0; } } - return matrix; + std::cout << "reached" << std::endl; } + std::vector> get() { return matrix_; } + private: - std::vector &)>> foos; - std::vector v; - Tape *tape; - unsigned int domain, codomain; + std::function(std::vector &)> foo_; + std::vector v_; + Tape *tape_; + unsigned int domain_, codomain_; + std::vector> matrix_; }; } // namespace xad diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 3dc54ed..fed388f 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -59,32 +59,76 @@ T foo(std::vector &x) return c + d + e; } +template +T single(std::vector &x) +{ + return x[0] * x[0] * x[0] + x[0]; +} + TEST(HessianTest, QuadraticForwardAdjoint) { - typedef xad::AReal> AD; + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; std::vector x = {3, 2}; - xad::Hessian h(quad, x); + xad::Hessian hes(quad, x, &tape); std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> computed_hessian = h.compute(); + std::vector> computed_hessian = hes.get(); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); } -TEST(HessianTest, ThreeVarQuadratic) +TEST(HessianTest, SingleInputForwardAdjoint) { typedef xad::fwd_adj mode; typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {3}; + xad::Hessian hes(single, x, &tape); + + std::vector> cross_hessian = {{18.0}}; + std::vector> computed_hessian = hes.get(); + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, QuadraticForwardForward) +{ + typedef xad::fwd_fwd mode; + typedef mode::active_type AD; + + std::vector x = {3, 2}; + xad::Hessian hes(quad, x); + + std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::vector> computed_hessian = hes.get(); + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, SingleInputForwardForward) +{ + typedef xad::fwd_fwd mode; + typedef mode::active_type AD; - std::vector x = {3, 2, 4}; - xad::Hessian hes(tquad, x); + std::vector x = {3}; + xad::Hessian hes(single, x); - std::vector> cross_hessian = { - {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}}; - std::vector> computed_hessian = hes.compute(); + std::vector> cross_hessian = {{18.0}}; + std::vector> computed_hessian = hes.get(); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index ad75d06..bacc040 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -41,6 +41,12 @@ T foo2(std::vector &x) return x[1] + sin(x[0]); } +template +std::vector foo(std::vector &x) +{ + return {foo1(x), foo2(x)}; +} + TEST(JacobianTest, SimpleJacobianAdjoint) { typedef xad::adj mode; @@ -49,15 +55,11 @@ TEST(JacobianTest, SimpleJacobianAdjoint) tape_type tape; - std::vector &)>> funcs(2); - funcs[0] = foo1; - funcs[1] = foo2; - std::vector x = {-2, 1}; - xad::Jacobian jac(funcs, x, &tape); + xad::Jacobian jac(foo, x, &tape); std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; - std::vector> computed_jacobian = jac.compute(); + std::vector> computed_jacobian = jac.get(); for (unsigned int i = 0; i < cross_jacobian.size(); i++) for (unsigned int j = 0; j < cross_jacobian.size(); j++) From e48fcf731bffe3e4776a1068cbeb1c840ed85e07 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Fri, 7 Jun 2024 11:04:19 +0100 Subject: [PATCH 07/30] Hessian & Jacobian now take matrix iterator (& corresponding tests) --- src/XAD/Hessian.hpp | 97 +++++++++++++++++++++++++-- src/XAD/Jacobian.hpp | 144 +++++++++++++++++++++++++++++++++++------ test/Hessian_test.cpp | 35 ++++++++++ test/Jacobian_test.cpp | 20 +++++- 4 files changed, 271 insertions(+), 25 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 6918941..f7f923a 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -32,7 +32,7 @@ template class Hessian { public: - // fwd_adj constructor + // fwd_adj 2d vector constructor Hessian(std::function &)> func, const std::vector &v, xad::Tape> *tape) : foo_(func), @@ -43,7 +43,20 @@ class Hessian compute(tape); } - // tapeless (fwd_fwd) constructor + // fwd_adj iterator constructor + template + Hessian(std::function &)> func, const std::vector &v, + xad::Tape> *tape, RowIterator first, RowIterator last) + : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) + { + if (std::distance(first, last) != domain_) + { + // throw exception + } + compute(tape, first, last); + } + + // tapeless (fwd_fwd) 2d vector constructor Hessian(std::function &)> func, std::vector &v) : foo_(func), v_(v), @@ -53,7 +66,20 @@ class Hessian compute(); } - // fwd_adj + // tapeless (fwd_fwd) iterator constructor + template + Hessian(std::function &)> func, std::vector &v, RowIterator first, + RowIterator last) + : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) + { + if (std::distance(first, last) != domain_) + { + // throw exception + } + compute(first, last); + } + + // fwd_adj 2d vector void compute(xad::Tape> *tape) { tape->registerInputs(v_); @@ -80,7 +106,39 @@ class Hessian } } - // fwd_fwd + // fwd_adj iterator + template + void compute(xad::Tape> *tape, RowIterator first, RowIterator last) + { + tape->registerInputs(v_); + + auto row = first; + + for (unsigned int i = 0; i < domain_; i++, row++) + { + derivative(value(v_[i])) = 1.0; + tape->newRecording(); + + T y = foo_(v_); + tape->registerOutput(y); + value(derivative(y)) = 1.0; + + tape->computeAdjoints(); + + auto col = row->begin(); + + for (unsigned int j = 0; j < domain_; j++, col++) + { + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) + // << "\n"; + *col = derivative(derivative(v_[j])); + } + + derivative(value(v_[i])) = 0.0; + } + } + + // fwd_fwd 2d vector void compute() { for (unsigned int i = 0; i < domain_; i++) @@ -94,7 +152,7 @@ class Hessian T y = foo_(v_); // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) - // << "\n"; + // << "\n";ß matrix_[i][j] = derivative(derivative(y)); derivative(value(v_[j])) = 0.0; @@ -104,6 +162,35 @@ class Hessian } } + // fwd_fwd iterator + template + void compute(RowIterator first, RowIterator last) + { + auto row = first; + + for (unsigned int i = 0; i < domain_; i++, row++) + { + value(derivative(v_[i])) = 1.0; + + auto col = row->begin(); + + for (unsigned int j = 0; j < domain_; j++, col++) + { + derivative(value(v_[j])) = 1.0; + + T y = foo_(v_); + + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) + // << "\n"; + *col = derivative(derivative(y)); + + derivative(value(v_[j])) = 0.0; + } + + value(derivative(v_[i])) = 0.0; + } + } + std::vector> get() { return matrix_; } private: diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index fffed58..0df9cee 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -32,45 +32,154 @@ template class Jacobian { public: + // adj 2d vector constructor Jacobian(std::function(std::vector &)> foo, const std::vector &v, xad::Tape *tape) - : foo_(foo), - v_(v), - tape_(tape), - domain_(static_cast(v_.size())), - codomain_(0), - matrix_(0) + : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) + { + compute(tape); + } + + // adj iterator constructor + template + Jacobian(std::function(std::vector &)> foo, const std::vector &v, + xad::Tape *tape, RowIterator first, RowIterator last) + : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) + { + if (std::distance(first, last) != domain_) + { + // throw exception + } + compute(tape, first, last); + } + + // tapeless fwd 2d vector constructor + Jacobian(std::function(std::vector &)> foo, const std::vector &v) + : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) { compute(); } + // tapeless fwd iterator constructor + template + Jacobian(std::function(std::vector &)> foo, const std::vector &v, + RowIterator first, RowIterator last) + : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) + { + if (std::distance(first, last) != domain_) + { + // throw exception + } + compute(first, last); + } + + // adj 2d vector + void compute(xad::Tape *tape) + { + tape->registerInputs(v_); + tape->newRecording(); + std::vector res = foo_(v_); + codomain_ = static_cast(res.size()); + matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); + + for (unsigned int i = 0; i < codomain_; i++) + { + tape->registerOutput(res[i]); + derivative(res[i]) = 1.0; + + tape->computeAdjoints(); + + for (unsigned int j = 0; j < domain_; j++) + { + std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[j]) << std::endl; + matrix_[i][j] = derivative(v_[j]); + } + + derivative(res[i]) = 0.0; + } + } + + // adj iterator + template + void compute(xad::Tape *tape, RowIterator first, RowIterator last) + { + tape->registerInputs(v_); + tape->newRecording(); + std::vector res = foo_(v_); + codomain_ = static_cast(res.size()); + auto row = first; + + for (unsigned int i = 0; i < codomain_; i++, row++) + { + tape->registerOutput(res[i]); + derivative(res[i]) = 1.0; + + tape->computeAdjoints(); + + auto col = row->begin(); + + for (unsigned int j = 0; j < domain_; j++, col++) + { + std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[j]) << std::endl; + *col = derivative(v_[j]); + } + + derivative(res[i]) = 0.0; + } + } + + // fwd 2d vector void compute() { - tape_->registerInputs(v_); std::vector res = foo_(v_); codomain_ = static_cast(res.size()); matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); for (unsigned int i = 0; i < domain_; i++) { + derivative(v_[i]) = 1.0; + for (unsigned int j = 0; j < codomain_; j++) { - derivative(v_[i]) = 1.0; - tape_->newRecording(); - T y = res[j]; - tape_->registerOutput(y); - derivative(y) = 1.0; - tape_->computeAdjoints(); + std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; - matrix_[i][j] = derivative(v_[i]); - derivative(v_[i]) = 0.0; + matrix_[i][j] = derivative(y); + // derivative(y) = 0.0; } + + derivative(v_[i]) = 0.0; } + } - std::cout << "reached" << std::endl; + // fwd iterator + template + void compute(RowIterator first, RowIterator last) + { + std::vector res = foo_(v_); + codomain_ = static_cast(res.size()); + + auto row = first; + + for (unsigned int i = 0; i < domain_; i++, row++) + { + derivative(v_[i]) = 1.0; + + auto col = row->begin(); + + for (unsigned int j = 0; j < codomain_; j++, col++) + { + T y = res[j]; + + std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; + + *col = derivative(y); + // derivative(y) = 0.0; + } + + derivative(v_[i]) = 0.0; + } } std::vector> get() { return matrix_; } @@ -78,7 +187,6 @@ class Jacobian private: std::function(std::vector &)> foo_; std::vector v_; - Tape *tape_; unsigned int domain_, codomain_; std::vector> matrix_; }; diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index fed388f..3ad9c86 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -84,6 +84,25 @@ TEST(HessianTest, QuadraticForwardAdjoint) ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); } +TEST(HessianTest, QuadraticForwardAdjointWithIterator) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {3, 2}; + std::vector> computed_hessian(x.size(), std::vector(x.size(), 0.0)); + xad::Hessian hes(quad, x, &tape, computed_hessian.begin(), computed_hessian.end()); + + std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} + TEST(HessianTest, SingleInputForwardAdjoint) { typedef xad::fwd_adj mode; @@ -119,6 +138,22 @@ TEST(HessianTest, QuadraticForwardForward) ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); } +TEST(HessianTest, QuadraticForwardForwardWithIterator) +{ + typedef xad::fwd_fwd mode; + typedef mode::active_type AD; + + std::vector x = {3, 2}; + std::vector> computed_hessian(x.size(), std::vector(x.size(), 0.0)); + xad::Hessian hes(quad, x, computed_hessian.begin(), computed_hessian.end()); + + std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + for (unsigned int i = 0; i < cross_hessian.size(); i++) + for (unsigned int j = 0; j < cross_hessian.size(); j++) + ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); +} + TEST(HessianTest, SingleInputForwardForward) { typedef xad::fwd_fwd mode; diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index bacc040..a0feaaa 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -47,7 +47,7 @@ std::vector foo(std::vector &x) return {foo1(x), foo2(x)}; } -TEST(JacobianTest, SimpleJacobianAdjoint) +TEST(JacobianTest, SimpleAdjoint) { typedef xad::adj mode; typedef mode::tape_type tape_type; @@ -55,9 +55,25 @@ TEST(JacobianTest, SimpleJacobianAdjoint) tape_type tape; - std::vector x = {-2, 1}; + std::vector x = {2, 1}; xad::Jacobian jac(foo, x, &tape); + std::vector> cross_jacobian = {{1.0, cos(2)}, {cos(1), 1.0}}; + std::vector> computed_jacobian = jac.get(); + + for (unsigned int i = 0; i < cross_jacobian.size(); i++) + for (unsigned int j = 0; j < cross_jacobian.size(); j++) + ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); +} + +TEST(JacobianTest, SimpleForward) +{ + typedef xad::fwd mode; + typedef mode::active_type AD; + + std::vector x = {-2, 1}; + xad::Jacobian jac(foo, x); + std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; std::vector> computed_jacobian = jac.get(); From a47b248ed213f638ca820e4a512bb1c25bef2f28 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Fri, 7 Jun 2024 11:42:52 +0100 Subject: [PATCH 08/30] Jacobian fwd mode & corresponding tests --- src/XAD/Jacobian.hpp | 21 +++++++++------------ test/Jacobian_test.cpp | 39 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 46 insertions(+), 14 deletions(-) diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 0df9cee..f9a3736 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -84,8 +84,9 @@ class Jacobian for (unsigned int i = 0; i < codomain_; i++) { - tape->registerOutput(res[i]); - derivative(res[i]) = 1.0; + T y = res[i]; + tape->registerOutput(y); + derivative(y) = 1.0; tape->computeAdjoints(); @@ -94,8 +95,7 @@ class Jacobian std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[j]) << std::endl; matrix_[i][j] = derivative(v_[j]); } - - derivative(res[i]) = 0.0; + derivative(y) = 0.0; } } @@ -131,8 +131,8 @@ class Jacobian // fwd 2d vector void compute() { - std::vector res = foo_(v_); - codomain_ = static_cast(res.size()); + // TODO: optimise to use less computations of foo_() + codomain_ = static_cast(foo_(v_).size()); matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); for (unsigned int i = 0; i < domain_; i++) @@ -141,12 +141,11 @@ class Jacobian for (unsigned int j = 0; j < codomain_; j++) { - T y = res[j]; + T y = foo_(v_)[j]; std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; matrix_[i][j] = derivative(y); - // derivative(y) = 0.0; } derivative(v_[i]) = 0.0; @@ -157,8 +156,7 @@ class Jacobian template void compute(RowIterator first, RowIterator last) { - std::vector res = foo_(v_); - codomain_ = static_cast(res.size()); + codomain_ = static_cast(foo_(v_).size()); auto row = first; @@ -170,12 +168,11 @@ class Jacobian for (unsigned int j = 0; j < codomain_; j++, col++) { - T y = res[j]; + T y = foo_(v_)[j]; std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; *col = derivative(y); - // derivative(y) = 0.0; } derivative(v_[i]) = 0.0; diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index a0feaaa..ee4bce0 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -55,10 +55,10 @@ TEST(JacobianTest, SimpleAdjoint) tape_type tape; - std::vector x = {2, 1}; + std::vector x = {3, 1}; xad::Jacobian jac(foo, x, &tape); - std::vector> cross_jacobian = {{1.0, cos(2)}, {cos(1), 1.0}}; + std::vector> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; std::vector> computed_jacobian = jac.get(); for (unsigned int i = 0; i < cross_jacobian.size(); i++) @@ -66,6 +66,25 @@ TEST(JacobianTest, SimpleAdjoint) ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); } +TEST(JacobianTest, SimpleAdjointIterator) +{ + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; + + std::vector x = {3, 1}; + std::vector> computed_jacobian(2.0, std::vector(2.0, 0.0)); + xad::Jacobian jac(foo, x, &tape, computed_jacobian.begin(), computed_jacobian.end()); + + std::vector> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; + + for (unsigned int i = 0; i < cross_jacobian.size(); i++) + for (unsigned int j = 0; j < cross_jacobian.size(); j++) + ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); +} + TEST(JacobianTest, SimpleForward) { typedef xad::fwd mode; @@ -81,3 +100,19 @@ TEST(JacobianTest, SimpleForward) for (unsigned int j = 0; j < cross_jacobian.size(); j++) ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); } + +TEST(JacobianTest, SimpleForwardIterator) +{ + typedef xad::fwd mode; + typedef mode::active_type AD; + + std::vector x = {-2, 1}; + std::vector> computed_jacobian(2.0, std::vector(2.0, 0.0)); + xad::Jacobian jac(foo, x, computed_jacobian.begin(), computed_jacobian.end()); + + std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; + + for (unsigned int i = 0; i < cross_jacobian.size(); i++) + for (unsigned int j = 0; j < cross_jacobian.size(); j++) + ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); +} From d619a9b24099e41467172963cd9783f9baabd0b7 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Fri, 7 Jun 2024 12:23:05 +0100 Subject: [PATCH 09/30] Jacobian adj mode implemented --- src/XAD/Jacobian.hpp | 61 ++++++++++++++++++-------------------------- 1 file changed, 25 insertions(+), 36 deletions(-) diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index f9a3736..78a9f18 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -27,7 +27,7 @@ namespace xad { - +// TODO: optimise to use less computations of foo_() for all compute() methods template class Jacobian { @@ -77,25 +77,22 @@ class Jacobian void compute(xad::Tape *tape) { tape->registerInputs(v_); - tape->newRecording(); - std::vector res = foo_(v_); - codomain_ = static_cast(res.size()); + codomain_ = static_cast(foo_(v_).size()); matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); - for (unsigned int i = 0; i < codomain_; i++) + for (unsigned int i = 0; i < domain_; i++) { - T y = res[i]; - tape->registerOutput(y); - derivative(y) = 1.0; - - tape->computeAdjoints(); - - for (unsigned int j = 0; j < domain_; j++) + for (unsigned int j = 0; j < codomain_; j++) { - std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[j]) << std::endl; - matrix_[i][j] = derivative(v_[j]); + tape->newRecording(); + T y = foo_(v_)[j]; + tape->registerOutput(y); + derivative(y) = 1.0; + tape->computeAdjoints(); + + // std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; + matrix_[i][j] = derivative(v_[i]); } - derivative(y) = 0.0; } } @@ -104,34 +101,30 @@ class Jacobian void compute(xad::Tape *tape, RowIterator first, RowIterator last) { tape->registerInputs(v_); - tape->newRecording(); - std::vector res = foo_(v_); - codomain_ = static_cast(res.size()); + codomain_ = static_cast(foo_(v_).size()); auto row = first; - for (unsigned int i = 0; i < codomain_; i++, row++) + for (unsigned int i = 0; i < domain_; i++, row++) { - tape->registerOutput(res[i]); - derivative(res[i]) = 1.0; - - tape->computeAdjoints(); - auto col = row->begin(); - for (unsigned int j = 0; j < domain_; j++, col++) + for (unsigned int j = 0; j < codomain_; j++, col++) { - std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[j]) << std::endl; - *col = derivative(v_[j]); - } + tape->newRecording(); + T y = foo_(v_)[j]; + tape->registerOutput(y); + derivative(y) = 1.0; + tape->computeAdjoints(); - derivative(res[i]) = 0.0; + // std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; + *col = derivative(v_[i]); + } } } // fwd 2d vector void compute() { - // TODO: optimise to use less computations of foo_() codomain_ = static_cast(foo_(v_).size()); matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); @@ -142,9 +135,7 @@ class Jacobian for (unsigned int j = 0; j < codomain_; j++) { T y = foo_(v_)[j]; - - std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - + // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; matrix_[i][j] = derivative(y); } @@ -169,9 +160,7 @@ class Jacobian for (unsigned int j = 0; j < codomain_; j++, col++) { T y = foo_(v_)[j]; - - std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - + // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; *col = derivative(y); } From 303ed3fceaf3159c8b78ec9d69fa669c191362c5 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Fri, 7 Jun 2024 12:44:18 +0100 Subject: [PATCH 10/30] Exception handling for RowIterator overloads & doc files added --- docs/README.md | 2 ++ docs/ref/hessian.md | 1 + docs/ref/jacobian.md | 1 + src/XAD/Hessian.hpp | 26 ++++++++++++++++++++------ src/XAD/Jacobian.hpp | 26 ++++++++++++++++++++------ 5 files changed, 44 insertions(+), 12 deletions(-) create mode 100644 docs/ref/hessian.md create mode 100644 docs/ref/jacobian.md diff --git a/docs/README.md b/docs/README.md index 02d334e..456ae70 100644 --- a/docs/README.md +++ b/docs/README.md @@ -15,4 +15,6 @@ which is imported to the website at [auto-differentiation.github.io](https://aut * [Tape](ref/tape.md) * [CheckpointCallback](ref/chkpt_cb.md) * [Exceptions](ref/exceptions.md) +* [Hessian](ref/hessian.md) +* [Jacobian](ref/jacobian.md) * [Version Information](ref/version.md) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md new file mode 100644 index 0000000..3eacf7b --- /dev/null +++ b/docs/ref/hessian.md @@ -0,0 +1 @@ +# Hessian \ No newline at end of file diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md new file mode 100644 index 0000000..36b54ff --- /dev/null +++ b/docs/ref/jacobian.md @@ -0,0 +1 @@ +# Jacobian \ No newline at end of file diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index f7f923a..25958c3 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -23,6 +23,7 @@ ******************************************************************************/ #include +#include #include namespace xad @@ -50,9 +51,10 @@ class Hessian : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) { if (std::distance(first, last) != domain_) - { - // throw exception - } + throw OutOfRange("Iterator not allocated enough space"); + static_assert(has_begin::value, + "RowIterator must dereference to a type that implements a begin() method"); + compute(tape, first, last); } @@ -73,9 +75,10 @@ class Hessian : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) { if (std::distance(first, last) != domain_) - { - // throw exception - } + throw OutOfRange("Iterator not allocated enough space"); + static_assert(has_begin::value, + "RowIterator must dereference to a type that implements a begin() method"); + compute(first, last); } @@ -198,5 +201,16 @@ class Hessian std::vector v_; unsigned int domain_; std::vector> matrix_; + + template + static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); + + template + static std::false_type has_begin_impl(...); + + template + struct has_begin : decltype(has_begin_impl(0)) + { + }; }; } // namespace xad diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 78a9f18..24edf4e 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -23,6 +23,7 @@ ******************************************************************************/ #include +#include #include namespace xad @@ -47,9 +48,10 @@ class Jacobian : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) { if (std::distance(first, last) != domain_) - { - // throw exception - } + throw OutOfRange("Iterator not allocated enough space"); + static_assert(has_begin::value, + "RowIterator must dereference to a type that implements a begin() method"); + compute(tape, first, last); } @@ -67,9 +69,10 @@ class Jacobian : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) { if (std::distance(first, last) != domain_) - { - // throw exception - } + throw OutOfRange("Iterator not allocated enough space"); + static_assert(has_begin::value, + "RowIterator must dereference to a type that implements a begin() method"); + compute(first, last); } @@ -175,5 +178,16 @@ class Jacobian std::vector v_; unsigned int domain_, codomain_; std::vector> matrix_; + + template + static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); + + template + static std::false_type has_begin_impl(...); + + template + struct has_begin : decltype(has_begin_impl(0)) + { + }; }; } // namespace xad From d2acd8a9522b36cfbdbf5b689e04862e4f5b6014 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Fri, 7 Jun 2024 13:19:12 +0100 Subject: [PATCH 11/30] Changed iterator tests to use std::list & changed CI to use msvc 14.4 for Windows --- .github/workflows/ci.yml | 2 +- test/Hessian_test.cpp | 44 +++++++++++++++++++++++++++++----------- test/Jacobian_test.cpp | 44 +++++++++++++++++++++++++++++----------- 3 files changed, 65 insertions(+), 25 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 82f8406..0724d55 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -6,7 +6,7 @@ jobs: strategy: fail-fast: false matrix: - toolset: ["14.0", "14.1", "14.2", "14.3"] + toolset: ["14.0", "14.1", "14.2", "14.4"] compiler: [msvc, clang] config: [Release, Debug] include: diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 3ad9c86..c4841e9 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -93,14 +93,24 @@ TEST(HessianTest, QuadraticForwardAdjointWithIterator) tape_type tape; std::vector x = {3, 2}; - std::vector> computed_hessian(x.size(), std::vector(x.size(), 0.0)); + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); xad::Hessian hes(quad, x, &tape, computed_hessian.begin(), computed_hessian.end()); - std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + auto row1 = computed_hessian.begin(), row2 = cross_hessian.begin(); + while (row1 != computed_hessian.end() && row2 != cross_hessian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } } TEST(HessianTest, SingleInputForwardAdjoint) @@ -144,14 +154,24 @@ TEST(HessianTest, QuadraticForwardForwardWithIterator) typedef mode::active_type AD; std::vector x = {3, 2}; - std::vector> computed_hessian(x.size(), std::vector(x.size(), 0.0)); + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); xad::Hessian hes(quad, x, computed_hessian.begin(), computed_hessian.end()); - std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + auto row1 = computed_hessian.begin(), row2 = cross_hessian.begin(); + while (row1 != computed_hessian.end() && row2 != cross_hessian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } } TEST(HessianTest, SingleInputForwardForward) diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index ee4bce0..bbafc26 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -75,14 +75,24 @@ TEST(JacobianTest, SimpleAdjointIterator) tape_type tape; std::vector x = {3, 1}; - std::vector> computed_jacobian(2.0, std::vector(2.0, 0.0)); + std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); xad::Jacobian jac(foo, x, &tape, computed_jacobian.begin(), computed_jacobian.end()); - std::vector> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; - - for (unsigned int i = 0; i < cross_jacobian.size(); i++) - for (unsigned int j = 0; j < cross_jacobian.size(); j++) - ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); + std::list> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; + + auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } } TEST(JacobianTest, SimpleForward) @@ -107,12 +117,22 @@ TEST(JacobianTest, SimpleForwardIterator) typedef mode::active_type AD; std::vector x = {-2, 1}; - std::vector> computed_jacobian(2.0, std::vector(2.0, 0.0)); + std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); xad::Jacobian jac(foo, x, computed_jacobian.begin(), computed_jacobian.end()); - std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; - - for (unsigned int i = 0; i < cross_jacobian.size(); i++) - for (unsigned int j = 0; j < cross_jacobian.size(); j++) - ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); + std::list> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; + + auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } } From d7722ca1b9b376e1c43b84470fef4cc09f532299 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Sun, 9 Jun 2024 14:51:16 +0100 Subject: [PATCH 12/30] Hessian & Jacobian impl. from OO to functional --- docs/ref/hessian.md | 2 +- src/XAD/Hessian.hpp | 268 +++++++++++++++++++---------------------- src/XAD/Jacobian.hpp | 230 +++++++++++++++-------------------- test/Hessian_test.cpp | 20 +-- test/Jacobian_test.cpp | 13 +- 5 files changed, 242 insertions(+), 291 deletions(-) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 3eacf7b..1c0567f 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -1 +1 @@ -# Hessian \ No newline at end of file +# Hessian diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 25958c3..cbafbfa 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -29,188 +29,172 @@ namespace xad { -template -class Hessian +namespace detail { - public: - // fwd_adj 2d vector constructor - Hessian(std::function &)> func, const std::vector &v, - xad::Tape> *tape) - : foo_(func), - v_(v), - domain_(static_cast(v_.size())), - matrix_(domain_, std::vector(domain_, 0.0)) - { - compute(tape); - } +template +static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); - // fwd_adj iterator constructor - template - Hessian(std::function &)> func, const std::vector &v, - xad::Tape> *tape, RowIterator first, RowIterator last) - : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) - { - if (std::distance(first, last) != domain_) - throw OutOfRange("Iterator not allocated enough space"); - static_assert(has_begin::value, - "RowIterator must dereference to a type that implements a begin() method"); +template +static std::false_type has_begin_impl(...); - compute(tape, first, last); - } +template +struct has_begin : decltype(has_begin_impl(0)) +{ +}; +} // namespace detail + +// fwd_adj 2d vector +template +std::vector>>> computeHessian( + std::vector>> &v, + std::function>(std::vector>> &)> func, + xad::Tape> *tape) +{ + unsigned int domain(static_cast(v.size())); + std::vector>>> matrix( + domain, std::vector>>(domain, 0.0)); + std::function>(std::vector>> &)> foo = func; - // tapeless (fwd_fwd) 2d vector constructor - Hessian(std::function &)> func, std::vector &v) - : foo_(func), - v_(v), - domain_(static_cast(v_.size())), - matrix_(domain_, std::vector(domain_, 0.0)) - { - compute(); - } + tape->registerInputs(v); - // tapeless (fwd_fwd) iterator constructor - template - Hessian(std::function &)> func, std::vector &v, RowIterator first, - RowIterator last) - : foo_(func), v_(v), domain_(static_cast(v_.size())), matrix_(0) + for (unsigned int i = 0; i < domain; i++) { - if (std::distance(first, last) != domain_) - throw OutOfRange("Iterator not allocated enough space"); - static_assert(has_begin::value, - "RowIterator must dereference to a type that implements a begin() method"); + derivative(value(v[i])) = 1.0; + tape->newRecording(); - compute(first, last); - } + xad::AReal> y = foo(v); + tape->registerOutput(y); + value(derivative(y)) = 1.0; - // fwd_adj 2d vector - void compute(xad::Tape> *tape) - { - tape->registerInputs(v_); + tape->computeAdjoints(); - for (unsigned int i = 0; i < domain_; i++) + for (unsigned int j = 0; j < domain; j++) { - derivative(value(v_[i])) = 1.0; - tape->newRecording(); - - T y = foo_(v_); - tape->registerOutput(y); - value(derivative(y)) = 1.0; + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) + // << "\n"; + matrix[i][j] = derivative(derivative(v[j])); + } - tape->computeAdjoints(); + derivative(value(v[i])) = 0.0; + } - for (unsigned int j = 0; j < domain_; j++) - { - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) - // << "\n"; - matrix_[i][j] = derivative(derivative(v_[j])); - } + return matrix; +} - derivative(value(v_[i])) = 0.0; - } - } +// fwd_adj iterator +template +void computeHessian( + std::vector>> &v, + std::function>(std::vector>> &)> func, + xad::Tape> *tape, RowIterator first, RowIterator last) +{ + unsigned int domain(static_cast(v.size())); + std::function>(std::vector>> &)> foo = func; - // fwd_adj iterator - template - void compute(xad::Tape> *tape, RowIterator first, RowIterator last) - { - tape->registerInputs(v_); + if (std::distance(first, last) != domain) + throw OutOfRange("Iterator not allocated enough space"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); - auto row = first; + tape->registerInputs(v); - for (unsigned int i = 0; i < domain_; i++, row++) - { - derivative(value(v_[i])) = 1.0; - tape->newRecording(); + auto row = first; - T y = foo_(v_); - tape->registerOutput(y); - value(derivative(y)) = 1.0; + for (unsigned int i = 0; i < domain; i++, row++) + { + derivative(value(v[i])) = 1.0; + tape->newRecording(); - tape->computeAdjoints(); + xad::AReal> y = foo(v); + tape->registerOutput(y); + value(derivative(y)) = 1.0; - auto col = row->begin(); + tape->computeAdjoints(); - for (unsigned int j = 0; j < domain_; j++, col++) - { - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) - // << "\n"; - *col = derivative(derivative(v_[j])); - } + auto col = row->begin(); - derivative(value(v_[i])) = 0.0; + for (unsigned int j = 0; j < domain; j++, col++) + { + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) + // << "\n"; + *col = derivative(derivative(v[j])); } + + derivative(value(v[i])) = 0.0; } +} - // fwd_fwd 2d vector - void compute() - { - for (unsigned int i = 0; i < domain_; i++) - { - value(derivative(v_[i])) = 1.0; +// fwd_fwd 2d vector +template +std::vector>>> computeHessian( + std::vector>> &v, + std::function>(std::vector>> &)> func) +{ + unsigned int domain(static_cast(v.size())); + std::vector>>> matrix( + domain, std::vector>>(domain, 0.0)); + std::function>(std::vector>> &)> foo = func; - for (unsigned int j = 0; j < domain_; j++) - { - derivative(value(v_[j])) = 1.0; + for (unsigned int i = 0; i < domain; i++) + { + value(derivative(v[i])) = 1.0; - T y = foo_(v_); + for (unsigned int j = 0; j < domain; j++) + { + derivative(value(v[j])) = 1.0; - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) - // << "\n";ß - matrix_[i][j] = derivative(derivative(y)); + xad::FReal> y = foo(v); - derivative(value(v_[j])) = 0.0; - } + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) + // << "\n";ß + matrix[i][j] = derivative(derivative(y)); - value(derivative(v_[i])) = 0.0; + derivative(value(v[j])) = 0.0; } - } - // fwd_fwd iterator - template - void compute(RowIterator first, RowIterator last) - { - auto row = first; - - for (unsigned int i = 0; i < domain_; i++, row++) - { - value(derivative(v_[i])) = 1.0; + value(derivative(v[i])) = 0.0; + } - auto col = row->begin(); + return matrix; +} - for (unsigned int j = 0; j < domain_; j++, col++) - { - derivative(value(v_[j])) = 1.0; +// fwd_fwd iterator +template +void computeHessian( + std::vector>> &v, + std::function>(std::vector>> &)> func, + RowIterator first, RowIterator last) +{ + unsigned int domain(static_cast(v.size())); + std::function>(std::vector>> &)> foo = func; - T y = foo_(v_); + if (std::distance(first, last) != domain) + throw OutOfRange("Iterator not allocated enough space"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) - // << "\n"; - *col = derivative(derivative(y)); + auto row = first; - derivative(value(v_[j])) = 0.0; - } + for (unsigned int i = 0; i < domain; i++, row++) + { + value(derivative(v[i])) = 1.0; - value(derivative(v_[i])) = 0.0; - } - } + auto col = row->begin(); - std::vector> get() { return matrix_; } + for (unsigned int j = 0; j < domain; j++, col++) + { + derivative(value(v[j])) = 1.0; - private: - std::function &)> foo_; - std::vector v_; - unsigned int domain_; - std::vector> matrix_; + xad::FReal> y = foo(v); - template - static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); + // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) + // << "\n"; + *col = derivative(derivative(y)); - template - static std::false_type has_begin_impl(...); + derivative(value(v[j])) = 0.0; + } - template - struct has_begin : decltype(has_begin_impl(0)) - { - }; -}; + value(derivative(v[i])) = 0.0; + } +} } // namespace xad diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 24edf4e..673504a 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -28,166 +28,132 @@ namespace xad { -// TODO: optimise to use less computations of foo_() for all compute() methods -template -class Jacobian +// adj 2d vector +template +std::vector>> computeJacobian( + std::vector> &v, + std::function>(std::vector> &)> foo, xad::Tape *tape) { - public: - // adj 2d vector constructor - Jacobian(std::function(std::vector &)> foo, const std::vector &v, - xad::Tape *tape) - : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) - { - compute(tape); - } - - // adj iterator constructor - template - Jacobian(std::function(std::vector &)> foo, const std::vector &v, - xad::Tape *tape, RowIterator first, RowIterator last) - : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) - { - if (std::distance(first, last) != domain_) - throw OutOfRange("Iterator not allocated enough space"); - static_assert(has_begin::value, - "RowIterator must dereference to a type that implements a begin() method"); + tape->registerInputs(v); + unsigned int domain = static_cast(v.size()), + codomain = static_cast(foo(v).size()); - compute(tape, first, last); - } + std::vector>> matrix_( + std::vector>>(codomain, std::vector>(domain, 0.0))); - // tapeless fwd 2d vector constructor - Jacobian(std::function(std::vector &)> foo, const std::vector &v) - : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) + for (unsigned int i = 0; i < domain; i++) { - compute(); + for (unsigned int j = 0; j < codomain; j++) + { + tape->newRecording(); + xad::AReal y = foo(v)[j]; + tape->registerOutput(y); + derivative(y) = 1.0; + tape->computeAdjoints(); + + // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; + matrix_[i][j] = derivative(v[i]); + } } - // tapeless fwd iterator constructor - template - Jacobian(std::function(std::vector &)> foo, const std::vector &v, - RowIterator first, RowIterator last) - : foo_(foo), v_(v), domain_(static_cast(v_.size())), codomain_(0), matrix_(0) - { - if (std::distance(first, last) != domain_) - throw OutOfRange("Iterator not allocated enough space"); - static_assert(has_begin::value, - "RowIterator must dereference to a type that implements a begin() method"); + return matrix_; +} - compute(first, last); - } +// adj iterator +template +void computeJacobian(std::vector> &v, + std::function>(std::vector> &)> foo, + xad::Tape *tape, RowIterator first, RowIterator last) +{ + tape->registerInputs(v); + unsigned int domain = static_cast(v.size()), + codomain = static_cast(foo(v).size()); + + if (std::distance(first, last) != domain) + throw OutOfRange("Iterator not allocated enough space"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + + auto row = first; - // adj 2d vector - void compute(xad::Tape *tape) + for (unsigned int i = 0; i < domain; i++, row++) { - tape->registerInputs(v_); - codomain_ = static_cast(foo_(v_).size()); - matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); + auto col = row->begin(); - for (unsigned int i = 0; i < domain_; i++) + for (unsigned int j = 0; j < codomain; j++, col++) { - for (unsigned int j = 0; j < codomain_; j++) - { - tape->newRecording(); - T y = foo_(v_)[j]; - tape->registerOutput(y); - derivative(y) = 1.0; - tape->computeAdjoints(); - - // std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; - matrix_[i][j] = derivative(v_[i]); - } + tape->newRecording(); + xad::AReal y = foo(v)[j]; + tape->registerOutput(y); + derivative(y) = 1.0; + tape->computeAdjoints(); + + // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; + *col = derivative(v[i]); } } +} - // adj iterator - template - void compute(xad::Tape *tape, RowIterator first, RowIterator last) +// fwd 2d vector +template +std::vector>> computeJacobian( + std::vector> &v, + std::function>(std::vector> &)> foo) +{ + unsigned int domain = static_cast(v.size()), + codomain = static_cast(foo(v).size()); + std::vector>> matrix_( + std::vector>>(codomain, std::vector>(domain, 0.0))); + + for (unsigned int i = 0; i < domain; i++) { - tape->registerInputs(v_); - codomain_ = static_cast(foo_(v_).size()); - auto row = first; + derivative(v[i]) = 1.0; - for (unsigned int i = 0; i < domain_; i++, row++) + for (unsigned int j = 0; j < codomain; j++) { - auto col = row->begin(); - - for (unsigned int j = 0; j < codomain_; j++, col++) - { - tape->newRecording(); - T y = foo_(v_)[j]; - tape->registerOutput(y); - derivative(y) = 1.0; - tape->computeAdjoints(); - - // std::cout << "df" << j << "/dx" << i << " = " << derivative(v_[i]) << std::endl; - *col = derivative(v_[i]); - } + xad::FReal y = foo(v)[j]; + // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; + matrix_[i][j] = derivative(y); } + + derivative(v[i]) = 0.0; } - // fwd 2d vector - void compute() - { - codomain_ = static_cast(foo_(v_).size()); - matrix_ = std::vector>(codomain_, std::vector(domain_, 0.0)); + return matrix_; +} - for (unsigned int i = 0; i < domain_; i++) - { - derivative(v_[i]) = 1.0; +// fwd iterator +template +void computeJacobian(std::vector> &v, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last) +{ + unsigned int domain = static_cast(v.size()), + codomain = static_cast(foo(v).size()); + std::vector>> matrix_( + std::vector>>(codomain, std::vector>(domain, 0.0))); - for (unsigned int j = 0; j < codomain_; j++) - { - T y = foo_(v_)[j]; - // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - matrix_[i][j] = derivative(y); - } + if (std::distance(first, last) != domain) + throw OutOfRange("Iterator not allocated enough space"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); - derivative(v_[i]) = 0.0; - } - } + auto row = first; - // fwd iterator - template - void compute(RowIterator first, RowIterator last) + for (unsigned int i = 0; i < domain; i++, row++) { - codomain_ = static_cast(foo_(v_).size()); + derivative(v[i]) = 1.0; - auto row = first; + auto col = row->begin(); - for (unsigned int i = 0; i < domain_; i++, row++) + for (unsigned int j = 0; j < codomain; j++, col++) { - derivative(v_[i]) = 1.0; - - auto col = row->begin(); - - for (unsigned int j = 0; j < codomain_; j++, col++) - { - T y = foo_(v_)[j]; - // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - *col = derivative(y); - } - - derivative(v_[i]) = 0.0; + xad::FReal y = foo(v)[j]; + // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; + *col = derivative(y); } - } - - std::vector> get() { return matrix_; } - - private: - std::function(std::vector &)> foo_; - std::vector v_; - unsigned int domain_, codomain_; - std::vector> matrix_; - - template - static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); - template - static std::false_type has_begin_impl(...); - - template - struct has_begin : decltype(has_begin_impl(0)) - { - }; -}; + derivative(v[i]) = 0.0; + } +} } // namespace xad diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index c4841e9..106f61e 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -74,10 +75,9 @@ TEST(HessianTest, QuadraticForwardAdjoint) tape_type tape; std::vector x = {3, 2}; - xad::Hessian hes(quad, x, &tape); std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> computed_hessian = hes.get(); + auto computed_hessian = xad::computeHessian(x, quad, &tape); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) @@ -94,7 +94,8 @@ TEST(HessianTest, QuadraticForwardAdjointWithIterator) std::vector x = {3, 2}; std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); - xad::Hessian hes(quad, x, &tape, computed_hessian.begin(), computed_hessian.end()); + xad::computeHessian( + x, quad, &tape, begin(computed_hessian), end(computed_hessian)); std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; @@ -122,10 +123,10 @@ TEST(HessianTest, SingleInputForwardAdjoint) tape_type tape; std::vector x = {3}; - xad::Hessian hes(single, x, &tape); std::vector> cross_hessian = {{18.0}}; - std::vector> computed_hessian = hes.get(); + std::vector> computed_hessian = + xad::computeHessian(x, single, &tape); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) @@ -138,10 +139,9 @@ TEST(HessianTest, QuadraticForwardForward) typedef mode::active_type AD; std::vector x = {3, 2}; - xad::Hessian hes(quad, x); std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> computed_hessian = hes.get(); + std::vector> computed_hessian = xad::computeHessian(x, quad); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) @@ -155,7 +155,8 @@ TEST(HessianTest, QuadraticForwardForwardWithIterator) std::vector x = {3, 2}; std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); - xad::Hessian hes(quad, x, computed_hessian.begin(), computed_hessian.end()); + xad::computeHessian( + x, quad, begin(computed_hessian), end(computed_hessian)); std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; @@ -180,10 +181,9 @@ TEST(HessianTest, SingleInputForwardForward) typedef mode::active_type AD; std::vector x = {3}; - xad::Hessian hes(single, x); std::vector> cross_hessian = {{18.0}}; - std::vector> computed_hessian = hes.get(); + std::vector> computed_hessian = xad::computeHessian(x, single); for (unsigned int i = 0; i < cross_hessian.size(); i++) for (unsigned int j = 0; j < cross_hessian.size(); j++) diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index bbafc26..60ae159 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -56,10 +56,10 @@ TEST(JacobianTest, SimpleAdjoint) tape_type tape; std::vector x = {3, 1}; - xad::Jacobian jac(foo, x, &tape); std::vector> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; - std::vector> computed_jacobian = jac.get(); + std::vector> computed_jacobian = + xad::computeJacobian(x, foo, &tape); for (unsigned int i = 0; i < cross_jacobian.size(); i++) for (unsigned int j = 0; j < cross_jacobian.size(); j++) @@ -76,9 +76,10 @@ TEST(JacobianTest, SimpleAdjointIterator) std::vector x = {3, 1}; std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); - xad::Jacobian jac(foo, x, &tape, computed_jacobian.begin(), computed_jacobian.end()); std::list> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; + xad::computeJacobian( + x, foo, &tape, begin(computed_jacobian), end(computed_jacobian)); auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) @@ -101,10 +102,9 @@ TEST(JacobianTest, SimpleForward) typedef mode::active_type AD; std::vector x = {-2, 1}; - xad::Jacobian jac(foo, x); std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; - std::vector> computed_jacobian = jac.get(); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < cross_jacobian.size(); i++) for (unsigned int j = 0; j < cross_jacobian.size(); j++) @@ -118,9 +118,10 @@ TEST(JacobianTest, SimpleForwardIterator) std::vector x = {-2, 1}; std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); - xad::Jacobian jac(foo, x, computed_jacobian.begin(), computed_jacobian.end()); std::list> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; + xad::computeJacobian( + x, foo, begin(computed_jacobian), end(computed_jacobian)); auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) From 1d260d2a0df617a81f198ed3bd6570b3d4e1b11b Mon Sep 17 00:00:00 2001 From: imraneamri Date: Sun, 9 Jun 2024 15:01:27 +0100 Subject: [PATCH 13/30] Fixed naming convention --- src/XAD/Jacobian.hpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 673504a..a3a4bca 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -38,7 +38,7 @@ std::vector>> computeJacobian( unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); - std::vector>> matrix_( + std::vector>> matrix( std::vector>>(codomain, std::vector>(domain, 0.0))); for (unsigned int i = 0; i < domain; i++) @@ -52,11 +52,11 @@ std::vector>> computeJacobian( tape->computeAdjoints(); // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; - matrix_[i][j] = derivative(v[i]); + matrix[i][j] = derivative(v[i]); } } - return matrix_; + return matrix; } // adj iterator @@ -102,7 +102,7 @@ std::vector>> computeJacobian( { unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); - std::vector>> matrix_( + std::vector>> matrix( std::vector>>(codomain, std::vector>(domain, 0.0))); for (unsigned int i = 0; i < domain; i++) @@ -113,13 +113,13 @@ std::vector>> computeJacobian( { xad::FReal y = foo(v)[j]; // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - matrix_[i][j] = derivative(y); + matrix[i][j] = derivative(y); } derivative(v[i]) = 0.0; } - return matrix_; + return matrix; } // fwd iterator @@ -130,7 +130,7 @@ void computeJacobian(std::vector> &v, { unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); - std::vector>> matrix_( + std::vector>> matrix( std::vector>>(codomain, std::vector>(domain, 0.0))); if (std::distance(first, last) != domain) From 08668da8f2407120f35fcb4b45b191d795c3f7b1 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Sun, 9 Jun 2024 15:31:38 +0100 Subject: [PATCH 14/30] Removed redundant variable --- src/XAD/Hessian.hpp | 2 +- src/XAD/Jacobian.hpp | 4 +--- test/Hessian_test.cpp | 2 +- test/Jacobian_test.cpp | 2 +- 4 files changed, 4 insertions(+), 6 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index cbafbfa..5c48803 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -1,6 +1,6 @@ /******************************************************************************* - Implementation of Hessian computing methods. + Implementation of hessian computing methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index a3a4bca..3f7fa17 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -1,6 +1,6 @@ /******************************************************************************* - Implementation of Jacobian computing methods. + Implementation of jacobian computing methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. @@ -130,8 +130,6 @@ void computeJacobian(std::vector> &v, { unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); - std::vector>> matrix( - std::vector>>(codomain, std::vector>(domain, 0.0))); if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 106f61e..62689a6 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -1,6 +1,6 @@ /******************************************************************************* - Tests for xad::Hessian class. + Tests for hessian methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 60ae159..7c3f242 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -1,6 +1,6 @@ /******************************************************************************* - Tests for xad::Jacobian class. + Tests for jacobian computation methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. From dc344eee656738cf4df9692bcb15b8f1bdb75db1 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Mon, 10 Jun 2024 11:00:22 +0100 Subject: [PATCH 15/30] QOL changes and new tests - New Hessian, Jacobian & TypeTraits tests - Hessians & Jacobians no loner populated with xad values and defaults to underlying type - Tape auto sets if needed - Fixed usage of has_begin - Removed some duplicate code - Hessian & Jacobian are no longer default included in XAD.hpp --- src/XAD/Hessian.hpp | 122 +++++--------- src/XAD/Jacobian.hpp | 90 +++++------ src/XAD/TypeTraits.hpp | 45 ++++++ src/XAD/XAD.hpp | 2 - test/CMakeLists.txt | 1 + test/Hessian_test.cpp | 340 +++++++++++++++++++++++++++++++-------- test/Jacobian_test.cpp | 213 ++++++++++++++++++------ test/TypeTraits_test.cpp | 53 ++++++ 8 files changed, 610 insertions(+), 256 deletions(-) create mode 100644 src/XAD/TypeTraits.hpp create mode 100644 test/TypeTraits_test.cpp diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 5c48803..9aee4dd 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -22,6 +22,9 @@ ******************************************************************************/ +#include +#include + #include #include #include @@ -29,72 +32,43 @@ namespace xad { -namespace detail -{ -template -static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); - -template -static std::false_type has_begin_impl(...); - -template -struct has_begin : decltype(has_begin_impl(0)) -{ -}; -} // namespace detail - // fwd_adj 2d vector template -std::vector>>> computeHessian( - std::vector>> &v, - std::function>(std::vector>> &)> func, - xad::Tape> *tape) +std::vector> computeHessian( + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + xad::Tape> *tape = xad::Tape>::getActive()) { - unsigned int domain(static_cast(v.size())); - std::vector>>> matrix( - domain, std::vector>>(domain, 0.0)); - std::function>(std::vector>> &)> foo = func; - - tape->registerInputs(v); - - for (unsigned int i = 0; i < domain; i++) - { - derivative(value(v[i])) = 1.0; - tape->newRecording(); - - xad::AReal> y = foo(v); - tape->registerOutput(y); - value(derivative(y)) = 1.0; - - tape->computeAdjoints(); - - for (unsigned int j = 0; j < domain; j++) - { - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) - // << "\n"; - matrix[i][j] = derivative(derivative(v[j])); - } - - derivative(value(v[i])) = 0.0; - } - + std::vector> matrix(vec.size(), std::vector(vec.size(), 0.0)); + computeHessian(vec, foo, begin(matrix), end(matrix), tape); return matrix; } // fwd_adj iterator template void computeHessian( - std::vector>> &v, - std::function>(std::vector>> &)> func, - xad::Tape> *tape, RowIterator first, RowIterator last) + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + RowIterator first, RowIterator last, + xad::Tape> *tape = xad::Tape>::getActive()) { + auto v(vec); unsigned int domain(static_cast(v.size())); - std::function>(std::vector>> &)> foo = func; if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert(detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert( + !xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + + if (!tape) + { + std::unique_ptr>> t; + t = std::unique_ptr>>(new xad::Tape>()); + if (!t) + throw xad::NoTapeException(); + tape = t.get(); + } tape->registerInputs(v); @@ -126,52 +100,30 @@ void computeHessian( // fwd_fwd 2d vector template -std::vector>>> computeHessian( - std::vector>> &v, - std::function>(std::vector>> &)> func) +std::vector> computeHessian( + const std::vector>> &vec, + std::function>(std::vector>> &)> foo) { - unsigned int domain(static_cast(v.size())); - std::vector>>> matrix( - domain, std::vector>>(domain, 0.0)); - std::function>(std::vector>> &)> foo = func; - - for (unsigned int i = 0; i < domain; i++) - { - value(derivative(v[i])) = 1.0; - - for (unsigned int j = 0; j < domain; j++) - { - derivative(value(v[j])) = 1.0; - - xad::FReal> y = foo(v); - - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) - // << "\n";ß - matrix[i][j] = derivative(derivative(y)); - - derivative(value(v[j])) = 0.0; - } - - value(derivative(v[i])) = 0.0; - } - + std::vector> matrix(vec.size(), std::vector(vec.size(), 0.0)); + computeHessian(vec, foo, begin(matrix), end(matrix)); return matrix; } // fwd_fwd iterator template void computeHessian( - std::vector>> &v, - std::function>(std::vector>> &)> func, + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, RowIterator first, RowIterator last) { + auto v(vec); unsigned int domain(static_cast(v.size())); - std::function>(std::vector>> &)> foo = func; if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert(detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert( + !xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); auto row = first; diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 3f7fa17..e55eec1 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -22,6 +22,9 @@ ******************************************************************************/ +#include +#include + #include #include #include @@ -30,49 +33,44 @@ namespace xad { // adj 2d vector template -std::vector>> computeJacobian( - std::vector> &v, - std::function>(std::vector> &)> foo, xad::Tape *tape) +std::vector> computeJacobian( + const std::vector> &vec, + std::function>(std::vector> &)> foo, + xad::Tape *tape = xad::Tape::getActive()) { - tape->registerInputs(v); - unsigned int domain = static_cast(v.size()), - codomain = static_cast(foo(v).size()); - - std::vector>> matrix( - std::vector>>(codomain, std::vector>(domain, 0.0))); - - for (unsigned int i = 0; i < domain; i++) - { - for (unsigned int j = 0; j < codomain; j++) - { - tape->newRecording(); - xad::AReal y = foo(v)[j]; - tape->registerOutput(y); - derivative(y) = 1.0; - tape->computeAdjoints(); - - // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; - matrix[i][j] = derivative(v[i]); - } - } - + auto v(vec); + std::vector> matrix(foo(v).size(), std::vector(v.size(), 0.0)); + computeJacobian(vec, foo, begin(matrix), end(matrix), tape); return matrix; } // adj iterator template -void computeJacobian(std::vector> &v, +void computeJacobian(const std::vector> &vec, std::function>(std::vector> &)> foo, - xad::Tape *tape, RowIterator first, RowIterator last) + RowIterator first, RowIterator last, + xad::Tape *tape = xad::Tape::getActive()) { + auto v(vec); + + if (!tape) + { + std::unique_ptr> t; + t = std::unique_ptr>(new xad::Tape()); + if (!t) + throw xad::NoTapeException(); + tape = t.get(); + } + tape->registerInputs(v); unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert(detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert( + !xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); auto row = first; @@ -96,45 +94,31 @@ void computeJacobian(std::vector> &v, // fwd 2d vector template -std::vector>> computeJacobian( - std::vector> &v, +std::vector> computeJacobian( + const std::vector> &vec, std::function>(std::vector> &)> foo) { - unsigned int domain = static_cast(v.size()), - codomain = static_cast(foo(v).size()); - std::vector>> matrix( - std::vector>>(codomain, std::vector>(domain, 0.0))); - - for (unsigned int i = 0; i < domain; i++) - { - derivative(v[i]) = 1.0; - - for (unsigned int j = 0; j < codomain; j++) - { - xad::FReal y = foo(v)[j]; - // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; - matrix[i][j] = derivative(y); - } - - derivative(v[i]) = 0.0; - } - + auto v(vec); + std::vector> matrix(foo(v).size(), std::vector(v.size(), 0.0)); + computeJacobian(vec, foo, begin(matrix), end(matrix)); return matrix; } // fwd iterator template -void computeJacobian(std::vector> &v, +void computeJacobian(const std::vector> &vec, std::function>(std::vector> &)> foo, RowIterator first, RowIterator last) { + auto v(vec); unsigned int domain = static_cast(v.size()), codomain = static_cast(foo(v).size()); if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert(detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert( + !xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); auto row = first; diff --git a/src/XAD/TypeTraits.hpp b/src/XAD/TypeTraits.hpp new file mode 100644 index 0000000..c7628b3 --- /dev/null +++ b/src/XAD/TypeTraits.hpp @@ -0,0 +1,45 @@ +/******************************************************************************* + + Implementation of helper traits. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +namespace xad +{ + +namespace detail +{ + +template +static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::false_type{}); + +template +static std::true_type has_begin_impl(...); + +template +struct has_begin : decltype(has_begin_impl(0)) +{ +}; + +} // namespace detail +} // namespace xad diff --git a/src/XAD/XAD.hpp b/src/XAD/XAD.hpp index db912f0..0a1a310 100644 --- a/src/XAD/XAD.hpp +++ b/src/XAD/XAD.hpp @@ -35,9 +35,7 @@ #include #include #include -#include #include -#include #include #include #include diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0a03df3..8202c96 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -46,6 +46,7 @@ set(testfiles PartialRollback_test.cpp Hessian_test.cpp Jacobian_test.cpp + TypeTraits_test.cpp ) if (CMAKE_COMPILER_IS_GNUCC) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 62689a6..ea1ddce 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -22,50 +22,20 @@ ******************************************************************************/ +#include #include #include #include #include #include +#include #include #include #include using namespace ::testing; -template -T quad(std::vector &x) -{ - T c = x[0] * x[0]; - T d = x[1] * x[1]; - return c + d; -} - -template -T tquad(std::vector &x) -{ - T c = x[0] * x[0]; - T d = x[1] * x[1]; - T e = x[2] * x[2]; - return c + d + e; -} - -template -T foo(std::vector &x) -{ - T c = exp(x[0]); - T d = sin(x[1]); - T e = cos(x[2]); - return c + d + e; -} - -template -T single(std::vector &x) -{ - return x[0] * x[0] * x[0] + x[0]; -} - TEST(HessianTest, QuadraticForwardAdjoint) { typedef xad::fwd_adj mode; @@ -74,14 +44,18 @@ TEST(HessianTest, QuadraticForwardAdjoint) tape_type tape; - std::vector x = {3, 2}; + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, quad, &tape); + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } TEST(HessianTest, QuadraticForwardAdjointWithIterator) @@ -92,15 +66,19 @@ TEST(HessianTest, QuadraticForwardAdjointWithIterator) tape_type tape; - std::vector x = {3, 2}; - std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); - xad::computeHessian( - x, quad, &tape, begin(computed_hessian), end(computed_hessian)); + std::vector x = {3.0, 2.0}; - std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - auto row1 = computed_hessian.begin(), row2 = cross_hessian.begin(); - while (row1 != computed_hessian.end() && row2 != cross_hessian.end()) + std::list> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); + xad::computeHessian(x, foo, begin(computed_hessian), + end(computed_hessian), &tape); + + auto row1 = computed_hessian.begin(), row2 = expected_hessian.begin(); + while (row1 != computed_hessian.end() && row2 != expected_hessian.end()) { auto col1 = row1->begin(), col2 = row2->begin(); while (col1 != row1->end() && col2 != row2->end()) @@ -122,15 +100,18 @@ TEST(HessianTest, SingleInputForwardAdjoint) tape_type tape; - std::vector x = {3}; + std::vector x = {3.0}; + + // f(x) = x[0]^3 + x[0] + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] * x[0] + x[0]; }; - std::vector> cross_hessian = {{18.0}}; - std::vector> computed_hessian = - xad::computeHessian(x, single, &tape); + std::vector> expected_hessian = {{18.0}}; - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } TEST(HessianTest, QuadraticForwardForward) @@ -138,14 +119,19 @@ TEST(HessianTest, QuadraticForwardForward) typedef xad::fwd_fwd mode; typedef mode::active_type AD; - std::vector x = {3, 2}; + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> computed_hessian = xad::computeHessian(x, quad); + std::vector> expected_hessian = {{2.0, 0.0}, // + {0.0, 2.0}}; - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + auto computed_hessian = xad::computeHessian(x, foo); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } TEST(HessianTest, QuadraticForwardForwardWithIterator) @@ -153,15 +139,20 @@ TEST(HessianTest, QuadraticForwardForwardWithIterator) typedef xad::fwd_fwd mode; typedef mode::active_type AD; - std::vector x = {3, 2}; - std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); - xad::computeHessian( - x, quad, begin(computed_hessian), end(computed_hessian)); + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; + + std::list> expected_hessian = {{2.0, 0.0}, // + {0.0, 2.0}}; - std::list> cross_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); + xad::computeHessian(x, foo, begin(computed_hessian), + end(computed_hessian)); - auto row1 = computed_hessian.begin(), row2 = cross_hessian.begin(); - while (row1 != computed_hessian.end() && row2 != cross_hessian.end()) + auto row1 = computed_hessian.begin(), row2 = expected_hessian.begin(); + while (row1 != computed_hessian.end() && row2 != expected_hessian.end()) { auto col1 = row1->begin(), col2 = row2->begin(); while (col1 != row1->end() && col2 != row2->end()) @@ -180,12 +171,219 @@ TEST(HessianTest, SingleInputForwardForward) typedef xad::fwd_fwd mode; typedef mode::active_type AD; - std::vector x = {3}; + std::vector x = {3.0}; + + // f(x) = x[0]^3 + x[0] + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] * x[0] + x[0]; }; + + std::vector> expected_hessian = {{18.0}}; + + auto computed_hessian = xad::computeHessian(x, foo); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, QuadraticThreeVariablesForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0}; + + // f(x) = x[0]^2 + x[1]^2 + x[2]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1] + x[2] * x[2]; }; + + std::vector> expected_hessian = { + {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}}; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, ComplexFunctionForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0, 4.0}; + + // f(x) = x[0] * sin(x[1]) + x[2] * exp(x[3]) + auto foo = [](std::vector &x) -> AD { return x[0] * sin(x[1]) + x[2] * exp(x[3]); }; + + std::vector> expected_hessian = {{0.0, cos(x[1]), 0.0, 0.0}, + {cos(x[1]), -x[0] * sin(x[1]), 0.0, 0.0}, + {0.0, 0.0, 0.0, exp(x[3])}, + {0.0, 0.0, exp(x[3]), x[2] * exp(x[3])}}; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, FourthOrderPolynomialForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0}; + + // f(x) = x[0]^4 + x[1]^4 + x[2]^4 + auto foo = [](std::vector &x) -> AD { return pow(x[0], 4) + pow(x[1], 4) + pow(x[2], 4); }; + + std::vector> expected_hessian = {{12.0 * x[0] * x[0], 0.0, 0.0}, + {0.0, 12.0 * x[1] * x[1], 0.0}, + {0.0, 0.0, 12.0 * x[2] * x[2]}}; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, HigherOrderInteractionForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0}; + + // f(x) = x[0] * x[1] * x[2] + auto foo = [](std::vector &x) -> AD { return x[0] * x[1] * x[2]; }; + + std::vector> expected_hessian = {{0.0, x[2], x[1]}, // + {x[2], 0.0, x[0]}, + {x[1], x[0], 0.0}}; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, QuadraticFourVariablesForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0, 4.0}; + + // f(x) = x[0]^2 + x[1]^2 + x[2]^2 + x[3]^2 + auto foo = [](std::vector &x) -> AD + { return x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3]; }; + + std::vector> expected_hessian = {{2.0, 0.0, 0.0, 0.0}, // + {0.0, 2.0, 0.0, 0.0}, + {0.0, 0.0, 2.0, 0.0}, + {0.0, 0.0, 0.0, 2.0}}; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, LargeHessianForwardAdjoint) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x(16); + for (size_t i = 0; i < 16; ++i) x[i] = static_cast(i + 1); + + // f(x) = sum(x[i]^2) + sum(x[i] * x[j]), i < j + auto foo = [](std::vector &x) -> AD + { + AD result = 0.0; + for (size_t i = 0; i < x.size(); ++i) result += x[i] * x[i]; + for (size_t i = 0; i < x.size(); ++i) + for (size_t j = i + 1; j < x.size(); ++j) result += x[i] * x[j]; + return result; + }; + + std::vector> expected_hessian(16, std::vector(16, 1.0)); + for (size_t i = 0; i < 16; ++i) expected_hessian[i][i] = 2.0; + + auto computed_hessian = xad::computeHessian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, LargeHessianForwardForward) +{ + typedef xad::fwd_fwd mode; + typedef mode::active_type AD; + + std::vector x(16); + for (size_t i = 0; i < 16; ++i) x[i] = static_cast(i + 1); + + // f(x) = sum(x[i]^2) + sum(x[i] * x[j]), i < j + auto foo = [](std::vector &x) -> AD + { + AD result = 0.0; + for (size_t i = 0; i < x.size(); ++i) result += x[i] * x[i]; + for (size_t i = 0; i < x.size(); ++i) + for (size_t j = i + 1; j < x.size(); ++j) result += x[i] * x[j]; + return result; + }; + + std::vector> expected_hessian(16, std::vector(16, 1.0)); + for (size_t i = 0; i < 16; ++i) expected_hessian[i][i] = 2.0; + + auto computed_hessian = xad::computeHessian(x, foo); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, QuadraticForwardAdjointAutoTape) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; + + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::vector> cross_hessian = {{18.0}}; - std::vector> computed_hessian = xad::computeHessian(x, single); + auto computed_hessian = xad::computeHessian(x, foo); - for (unsigned int i = 0; i < cross_hessian.size(); i++) - for (unsigned int j = 0; j < cross_hessian.size(); j++) - ASSERT_EQ(cross_hessian[i][j], computed_hessian[i][j]); + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 7c3f242..6338cca 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -21,32 +21,17 @@ along with this program. If not, see . ******************************************************************************/ - +/* +#include #include #include #include #include +#include +#include #include -template -T foo1(std::vector &x) -{ - return x[0] + sin(x[1]); -} - -template -T foo2(std::vector &x) -{ - return x[1] + sin(x[0]); -} - -template -std::vector foo(std::vector &x) -{ - return {foo1(x), foo2(x)}; -} - TEST(JacobianTest, SimpleAdjoint) { typedef xad::adj mode; @@ -55,18 +40,22 @@ TEST(JacobianTest, SimpleAdjoint) tape_type tape; - std::vector x = {3, 1}; + std::vector x = {3.0, 1.0}; + + // f(x) = [ x[0] + sin(x[1]), x[1] + sin(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::vector> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; - std::vector> computed_jacobian = - xad::computeJacobian(x, foo, &tape); + std::vector> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; - for (unsigned int i = 0; i < cross_jacobian.size(); i++) - for (unsigned int j = 0; j < cross_jacobian.size(); j++) - ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); + auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); } -TEST(JacobianTest, SimpleAdjointIterator) +TEST(JacobianTest, SimpleAdjointIteratorAutoTape) { typedef xad::adj mode; typedef mode::tape_type tape_type; @@ -74,15 +63,20 @@ TEST(JacobianTest, SimpleAdjointIterator) tape_type tape; - std::vector x = {3, 1}; - std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); + std::vector x = {3.0, 1.0}; + + // f(x) = [ x[0] + sin(x[1]), x[1] + sin(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; + + std::list> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; - std::list> cross_jacobian = {{1.0, cos(x[0])}, {cos(x[1]), 1.0}}; + std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( - x, foo, &tape, begin(computed_jacobian), end(computed_jacobian)); + x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) + auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) { auto col1 = row1->begin(), col2 = row2->begin(); while (col1 != row1->end() && col2 != row2->end()) @@ -101,14 +95,19 @@ TEST(JacobianTest, SimpleForward) typedef xad::fwd mode; typedef mode::active_type AD; - std::vector x = {-2, 1}; + std::vector x = {-2.0, 1.0}; - std::vector> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; - std::vector> computed_jacobian = xad::computeJacobian(x, foo); + // f(x) = [ x[0] + sin(x[1]), x[1] + sin(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - for (unsigned int i = 0; i < cross_jacobian.size(); i++) - for (unsigned int j = 0; j < cross_jacobian.size(); j++) - ASSERT_EQ(cross_jacobian[i][j], computed_jacobian[i][j]); + std::vector> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + + auto computed_jacobian = xad::computeJacobian(x, foo); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); } TEST(JacobianTest, SimpleForwardIterator) @@ -116,15 +115,119 @@ TEST(JacobianTest, SimpleForwardIterator) typedef xad::fwd mode; typedef mode::active_type AD; - std::vector x = {-2, 1}; - std::list> computed_jacobian(2.0, std::list(2.0, 0.0)); + std::vector x = {-2.0, 1.0}; + + // f(x) = [ x[0] + sin(x[1]), x[1] + sin(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; + + std::list> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + + std::list> computed_jacobian(2, std::list(2, 0.0)); + xad::computeJacobian( + x, foo, begin(computed_jacobian), end(computed_jacobian)); + + auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } +} + +TEST(JacobianTest, ComplexFunctionAdjoint) +{ + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; + + std::vector x = {1.0, 2.0, 3.0, 4.0}; + + // f(x) = [ x[0] * x[1], x[2] * exp(x[3]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] * x[1], x[2] * exp(x[3])}; }; + + std::vector> expected_jacobian = {{x[1], x[0], 0.0, 0.0}, + {0.0, 0.0, exp(x[3]), x[2] * exp(x[3])}}; + + auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} + +TEST(JacobianTest, DomainLargerThanCodomainForward) +{ + typedef xad::fwd mode; + typedef mode::active_type AD; + + std::vector x = {1.0, 2.0, 3.0, 4.0}; + + // f(x) = [ x[0] + x[1], x[2] * x[3] ] + auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[2] * x[3]}; }; + + std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, {0.0, 0.0, x[3], x[2]}}; + + auto computed_jacobian = xad::computeJacobian(x, foo); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} + +TEST(JacobianTest, DomainSmallerThanCodomainAdjoint) +{ + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; + + std::vector x = {2.0, 3.0}; + + // f(x) = [ x[0] + x[1], x[0] - x[1], x[0] * x[1] ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + x[1], x[0] - x[1], x[0] * x[1]}; }; - std::list> cross_jacobian = {{1.0, cos(-2)}, {cos(1), 1.0}}; + std::vector> expected_jacobian = {{1.0, 1.0}, {1.0, -1.0}, {x[1], x[0]}}; + + auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} + +TEST(JacobianTest, ComplexDomainNotEqualCodomainForwardIterator) +{ + typedef xad::fwd mode; + typedef mode::active_type AD; + + std::vector x = {1.0, 2.0, 3.0}; + + // f(x) = [ x[0] + x[1], x[1] * x[2], exp(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + x[1], x[1] * x[2], exp(x[0])}; }; + + std::list> expected_jacobian = { + {1.0, 1.0, 0.0}, {0.0, x[2], x[1]}, {exp(x[0]), 0.0, 0.0}}; + + std::list> computed_jacobian(3, std::list(3, 0.0)); xad::computeJacobian( - x, foo, begin(computed_jacobian), end(computed_jacobian)); + x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = cross_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != cross_jacobian.end()) + auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) { auto col1 = row1->begin(), col2 = row2->begin(); while (col1 != row1->end() && col2 != row2->end()) @@ -137,3 +240,23 @@ TEST(JacobianTest, SimpleForwardIterator) row2++; } } + +TEST(JacobianTest, TrigonometricFunctionForward) +{ + typedef xad::fwd mode; + typedef mode::active_type AD; + + std::vector x = {M_PI / 4, M_PI / 3}; + + // f(x) = [ sin(x[0]), cos(x[1]) ] + auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; + + std::vector> expected_jacobian = {{cos(x[0]), 0.0}, {0.0, -sin(x[1])}}; + + auto computed_jacobian = xad::computeJacobian(x, foo); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} +*/ diff --git a/test/TypeTraits_test.cpp b/test/TypeTraits_test.cpp new file mode 100644 index 0000000..5b9dd20 --- /dev/null +++ b/test/TypeTraits_test.cpp @@ -0,0 +1,53 @@ +/******************************************************************************* + + Tests for helper traits. + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +#include +#include +#include + +TEST(TypeTraits, HasBeginVectorOfVectorsIterator) +{ + using It = std::vector>::iterator; + EXPECT_TRUE(xad::detail::has_begin::value); +} + +TEST(TypeTraits, HasBeginVectorOfVectors) +{ + using It = std::vector>; + EXPECT_FALSE(xad::detail::has_begin::value); +} + +TEST(TypeTraits, HasBeginListIterator) +{ + using It = std::list::iterator; + EXPECT_TRUE(xad::detail::has_begin::value); +} + +TEST(TypeTraits, HasBeginList) +{ + using It = std::list; + EXPECT_FALSE(xad::detail::has_begin::value); +} From 0a44e626e65789c6e5eb3e9d677c43ac641b33b9 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Mon, 10 Jun 2024 13:09:09 +0100 Subject: [PATCH 16/30] Test fixes, finished Jacobian impl., guard headers and has_begin reverted --- src/XAD/Hessian.hpp | 35 +++++++------------------- src/XAD/Jacobian.hpp | 53 ++++++++++++++++------------------------ src/XAD/TypeTraits.hpp | 5 ++-- test/Hessian_test.cpp | 3 --- test/Jacobian_test.cpp | 22 +++++++++++++++-- test/TypeTraits_test.cpp | 12 ++++++--- 6 files changed, 62 insertions(+), 68 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 9aee4dd..6fc7a34 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -22,6 +22,7 @@ ******************************************************************************/ +#pragma once #include #include @@ -52,48 +53,38 @@ void computeHessian( RowIterator first, RowIterator last, xad::Tape> *tape = xad::Tape>::getActive()) { - auto v(vec); - unsigned int domain(static_cast(v.size())); + unsigned int domain(static_cast(vec.size())); if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); static_assert( - !xad::detail::has_begin::value_type>::value, + xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); - + std::unique_ptr>> t; if (!tape) { - std::unique_ptr>> t; t = std::unique_ptr>>(new xad::Tape>()); - if (!t) - throw xad::NoTapeException(); tape = t.get(); } + auto v(vec); tape->registerInputs(v); auto row = first; - for (unsigned int i = 0; i < domain; i++, row++) { derivative(value(v[i])) = 1.0; tape->newRecording(); - xad::AReal> y = foo(v); tape->registerOutput(y); value(derivative(y)) = 1.0; - tape->computeAdjoints(); auto col = row->begin(); - for (unsigned int j = 0; j < domain; j++, col++) { - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(v[j])) - // << "\n"; *col = derivative(derivative(v[j])); } - derivative(value(v[i])) = 0.0; } } @@ -116,36 +107,28 @@ void computeHessian( std::function>(std::vector>> &)> foo, RowIterator first, RowIterator last) { - auto v(vec); - unsigned int domain(static_cast(v.size())); + unsigned int domain(static_cast(vec.size())); if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); static_assert( - !xad::detail::has_begin::value_type>::value, + xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); - auto row = first; + auto v(vec); + auto row = first; for (unsigned int i = 0; i < domain; i++, row++) { value(derivative(v[i])) = 1.0; - auto col = row->begin(); - for (unsigned int j = 0; j < domain; j++, col++) { derivative(value(v[j])) = 1.0; - xad::FReal> y = foo(v); - - // std::cout << "d2y/dx" << i << "dx" << j << " = " << derivative(derivative(y)) - // << "\n"; *col = derivative(derivative(y)); - derivative(value(v[j])) = 0.0; } - value(derivative(v[i])) = 0.0; } } diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index e55eec1..5dab2ce 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -22,6 +22,7 @@ ******************************************************************************/ +#pragma once #include #include @@ -52,42 +53,35 @@ void computeJacobian(const std::vector> &vec, xad::Tape *tape = xad::Tape::getActive()) { auto v(vec); + unsigned int domain = static_cast(vec.size()), + codomain = static_cast(foo(v).size()); + if (std::distance(first, last) != codomain) + throw OutOfRange("Iterator not allocated enough space"); + static_assert( + xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + std::unique_ptr> t; if (!tape) { - std::unique_ptr> t; t = std::unique_ptr>(new xad::Tape()); - if (!t) - throw xad::NoTapeException(); tape = t.get(); } tape->registerInputs(v); - unsigned int domain = static_cast(v.size()), - codomain = static_cast(foo(v).size()); - - if (std::distance(first, last) != domain) - throw OutOfRange("Iterator not allocated enough space"); - static_assert( - !xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); auto row = first; - - for (unsigned int i = 0; i < domain; i++, row++) + for (unsigned int i = 0; i < codomain; i++, row++) { auto col = row->begin(); - - for (unsigned int j = 0; j < codomain; j++, col++) + for (unsigned int j = 0; j < domain; j++, col++) { tape->newRecording(); - xad::AReal y = foo(v)[j]; + xad::AReal y = foo(v)[i]; tape->registerOutput(y); derivative(y) = 1.0; tape->computeAdjoints(); - - // std::cout << "df" << j << "/dx" << i << " = " << derivative(v[i]) << std::endl; - *col = derivative(v[i]); + *col = derivative(v[j]); } } } @@ -111,31 +105,26 @@ void computeJacobian(const std::vector> &vec, RowIterator first, RowIterator last) { auto v(vec); - unsigned int domain = static_cast(v.size()), + unsigned int domain = static_cast(vec.size()), codomain = static_cast(foo(v).size()); - if (std::distance(first, last) != domain) + if (std::distance(first, last) != codomain) throw OutOfRange("Iterator not allocated enough space"); static_assert( - !xad::detail::has_begin::value_type>::value, + xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); auto row = first; - - for (unsigned int i = 0; i < domain; i++, row++) + for (unsigned int i = 0; i < codomain; i++, row++) { - derivative(v[i]) = 1.0; - auto col = row->begin(); - - for (unsigned int j = 0; j < codomain; j++, col++) + for (unsigned int j = 0; j < domain; j++, col++) { - xad::FReal y = foo(v)[j]; - // std::cout << "df" << j << "/dx" << i << " = " << derivative(y) << std::endl; + derivative(v[j]) = 1.0; + xad::FReal y = foo(v)[i]; *col = derivative(y); + derivative(v[j]) = 0.0; } - - derivative(v[i]) = 0.0; } } } // namespace xad diff --git a/src/XAD/TypeTraits.hpp b/src/XAD/TypeTraits.hpp index c7628b3..a026a92 100644 --- a/src/XAD/TypeTraits.hpp +++ b/src/XAD/TypeTraits.hpp @@ -22,6 +22,7 @@ ******************************************************************************/ +#pragma once #include namespace xad @@ -31,10 +32,10 @@ namespace detail { template -static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::false_type{}); +static auto has_begin_impl(int) -> decltype(std::declval().begin(), std::true_type{}); template -static std::true_type has_begin_impl(...); +static std::false_type has_begin_impl(...); template struct has_begin : decltype(has_begin_impl(0)) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index ea1ddce..12aa2d8 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -370,9 +370,6 @@ TEST(HessianTest, QuadraticForwardAdjointAutoTape) { typedef xad::fwd_adj mode; typedef mode::active_type AD; - typedef mode::tape_type tape_type; - - tape_type tape; std::vector x = {3.0, 2.0}; diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 6338cca..2957394 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -21,7 +21,7 @@ along with this program. If not, see . ******************************************************************************/ -/* + #include #include @@ -259,4 +259,22 @@ TEST(JacobianTest, TrigonometricFunctionForward) for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); } -*/ + +TEST(JacobianTest, TrigonometricFunctionAdjointAutoTape) +{ + typedef xad::adj mode; + typedef mode::active_type AD; + + std::vector x = {M_PI / 4, M_PI / 3}; + + // f(x) = [ sin(x[0]), cos(x[1]) ] + auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; + + std::vector> expected_jacobian = {{cos(x[0]), 0.0}, {0.0, -sin(x[1])}}; + + auto computed_jacobian = xad::computeJacobian(x, foo); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} diff --git a/test/TypeTraits_test.cpp b/test/TypeTraits_test.cpp index 5b9dd20..c45af2d 100644 --- a/test/TypeTraits_test.cpp +++ b/test/TypeTraits_test.cpp @@ -31,23 +31,29 @@ TEST(TypeTraits, HasBeginVectorOfVectorsIterator) { using It = std::vector>::iterator; - EXPECT_TRUE(xad::detail::has_begin::value); + EXPECT_FALSE(xad::detail::has_begin::value); } TEST(TypeTraits, HasBeginVectorOfVectors) { using It = std::vector>; - EXPECT_FALSE(xad::detail::has_begin::value); + EXPECT_TRUE(xad::detail::has_begin::value); } TEST(TypeTraits, HasBeginListIterator) { using It = std::list::iterator; - EXPECT_TRUE(xad::detail::has_begin::value); + EXPECT_FALSE(xad::detail::has_begin::value); } TEST(TypeTraits, HasBeginList) { using It = std::list; + EXPECT_TRUE(xad::detail::has_begin::value); +} + +TEST(TypeTraits, HasBeginNonIterable) +{ + using It = size_t; EXPECT_FALSE(xad::detail::has_begin::value); } From 5c662071963c15ea7d2ae277ee40b365ea6469e4 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Mon, 10 Jun 2024 13:14:23 +0100 Subject: [PATCH 17/30] Added M_PI def for windows --- test/Jacobian_test.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 2957394..9bf6ddb 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -32,6 +32,10 @@ #include #include +#ifndef M_PI +#define M_PI 3.141592653589793238462643383279502884197169399 +#endif + TEST(JacobianTest, SimpleAdjoint) { typedef xad::adj mode; From a945864f3795fa42074ce25ed08794487f0ffed5 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Mon, 10 Jun 2024 14:16:07 +0100 Subject: [PATCH 18/30] Added MATH_DEFINES in tests for Windows & added basic documentation --- docs/ref/hessian.md | 21 +++++++++++++++++++++ docs/ref/jacobian.md | 23 ++++++++++++++++++++++- test/Jacobian_test.cpp | 5 +---- 3 files changed, 44 insertions(+), 5 deletions(-) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 1c0567f..b14c15b 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -1 +1,22 @@ # Hessian + +## Overview + +XAD implements a set of methods to compute the Hessian matrix of a function and its inputs in `XAD/Hessian.hpp`. + +Note that the Hessian header is not automatically included with `XAD/XAD.hpp`. +Users must include it as needed. + + +## Specialisations +Hessians can be computed in `fwd_adj` or `fwd_fwd` higher-order mode. + +The `computeHessian()` method takes a set of variables packaged in a `std::vector` and a function in the format `T foo(std::vector)`. + +If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. + +#### `fwd_adj` +This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. + +#### `fwd_fwd` +This mode does not require a Tape and can help reduce the overhead that comes with one. diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 36b54ff..14f27fa 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -1 +1,22 @@ -# Jacobian \ No newline at end of file +# Jacobian + +## Overview + +XAD implements a set of methods to compute the Jacobian matrix of a function and its inputs in `XAD/Jacobian.hpp`. + +Note that the Jacobian header is not automatically included with `XAD/XAD.hpp`. +Users must include it as needed. + + +## Specialisations +Jacobians can be computed in `adj` or `fwd` mode. + +The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function in the format `std::vector foo(std::vector)`. + +If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. + +#### `adj` +This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. + +#### `fwd_fwd` +This mode does not require a Tape and can help reduce the overhead that comes with one. diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 9bf6ddb..323fb3f 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -22,6 +22,7 @@ ******************************************************************************/ +#define _USE_MATH_DEFINES #include #include @@ -32,10 +33,6 @@ #include #include -#ifndef M_PI -#define M_PI 3.141592653589793238462643383279502884197169399 -#endif - TEST(JacobianTest, SimpleAdjoint) { typedef xad::adj mode; From d28cdfe715fd83a3fb86794b4aa4d4fbbc8aa817 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Mon, 10 Jun 2024 14:40:36 +0100 Subject: [PATCH 19/30] Pre-merge code cleanup --- docs/ref/hessian.md | 4 +++- docs/ref/jacobian.md | 4 +++- src/XAD/Hessian.hpp | 3 ++- src/XAD/Jacobian.hpp | 3 ++- src/XAD/TypeTraits.hpp | 2 +- 5 files changed, 11 insertions(+), 5 deletions(-) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index b14c15b..59ebdd7 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -9,14 +9,16 @@ Users must include it as needed. ## Specialisations + Hessians can be computed in `fwd_adj` or `fwd_fwd` higher-order mode. The `computeHessian()` method takes a set of variables packaged in a `std::vector` and a function in the format `T foo(std::vector)`. - If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. #### `fwd_adj` + This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` + This mode does not require a Tape and can help reduce the overhead that comes with one. diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 14f27fa..d443641 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -9,14 +9,16 @@ Users must include it as needed. ## Specialisations + Jacobians can be computed in `adj` or `fwd` mode. The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function in the format `std::vector foo(std::vector)`. - If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. #### `adj` + This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` + This mode does not require a Tape and can help reduce the overhead that comes with one. diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 6fc7a34..0667577 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -1,6 +1,6 @@ /******************************************************************************* - Implementation of hessian computing methods. + Implementation of Hessian matrix computing methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. @@ -23,6 +23,7 @@ ******************************************************************************/ #pragma once + #include #include diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 5dab2ce..bcaca9c 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -1,6 +1,6 @@ /******************************************************************************* - Implementation of jacobian computing methods. + Implementation of Jacobian matrix computing methods. This file is part of XAD, a comprehensive C++ library for automatic differentiation. @@ -23,6 +23,7 @@ ******************************************************************************/ #pragma once + #include #include diff --git a/src/XAD/TypeTraits.hpp b/src/XAD/TypeTraits.hpp index a026a92..17c2928 100644 --- a/src/XAD/TypeTraits.hpp +++ b/src/XAD/TypeTraits.hpp @@ -23,11 +23,11 @@ ******************************************************************************/ #pragma once + #include namespace xad { - namespace detail { From 93f5d90ef318301f15c4bde10ced2ae8755b99ff Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 14:20:33 +0100 Subject: [PATCH 20/30] Tests, doc & semantics: - Added sample programs for Hessian and Jacobian - Added doc for Hessian and Jacobian - Tweaked TypeTraits tests to be easier to understand - Changed order of some checks for iterators Upcoming: - More efficient fwd Jacobian computation - Possibly less evaluations of foo for Jacobian --- docs/ref/headers.md | 8 ++- docs/ref/hessian.md | 86 ++++++++++++++++++++++++++++++--- docs/ref/jacobian.md | 78 +++++++++++++++++++++++++++--- samples/CMakeLists.txt | 2 + samples/Hessian/CMakeLists.txt | 25 ++++++++++ samples/Hessian/main.cpp | 59 ++++++++++++++++++++++ samples/Jacobian/CMakeLists.txt | 25 ++++++++++ samples/Jacobian/main.cpp | 58 ++++++++++++++++++++++ src/XAD/Hessian.hpp | 4 ++ src/XAD/Jacobian.hpp | 47 ++++++++++-------- test/TypeTraits_test.cpp | 9 ++-- 11 files changed, 362 insertions(+), 39 deletions(-) create mode 100644 samples/Hessian/CMakeLists.txt create mode 100644 samples/Hessian/main.cpp create mode 100644 samples/Jacobian/CMakeLists.txt create mode 100644 samples/Jacobian/main.cpp diff --git a/docs/ref/headers.md b/docs/ref/headers.md index 63746d3..112968b 100644 --- a/docs/ref/headers.md +++ b/docs/ref/headers.md @@ -7,12 +7,16 @@ XAD provides a general header `XAD/XAD.hpp`, which includes all headers that are commonly needed to work with XAD. Typically, this is all that clients need to include. -There are two additional headers provided that can be included on demand: +There are four additional headers provided that can be included on demand: -* `XAD/Complex.hpp` - for using complex numbers with XAD data types (see [Complex](complex.md)). +* `XAD/Complex.hpp` - For using complex numbers with XAD data types (see [Complex](complex.md)). This header should be included wherever [`#!c++ std::complex`](https://en.cppreference.com/w/cpp/numeric/complex) is used. * `XAD/StdCompatibility.hpp` - This header imports the XAD math functions into the `std` namespace, for compatibility reasons. It enables using constructs like [`#!c++ std::sin(x)`](https://en.cppreference.com/w/cpp/numeric/math/sin) where `x` is an XAD type. Additionally, it also specialises [`#!c++ std::numeric_limits`](https://en.cppreference.com/w/cpp/types/numeric_limits) for the XAD data types, so that it provides traits similar to the standard floating point types. +* `XAD/Hessian.hpp` - Imports methods for computing the Hessian matrix of a single output function into the `xad` + namespace. +* `XAD/Jacobian.hpp` - Imports methods for computing the Jacobian matrix of a function with multiple inputs and multiple + outputs into the `xad` namespace. diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 59ebdd7..89761d8 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -2,23 +2,97 @@ ## Overview -XAD implements a set of methods to compute the Hessian matrix of a function and its inputs in `XAD/Hessian.hpp`. +XAD implements a set of methods to compute the Hessian matrix of a function in `XAD/Hessian.hpp`. Note that the Hessian header is not automatically included with `XAD/XAD.hpp`. Users must include it as needed. +Hessians can be computed in `fwd_adj` or `fwd_fwd` higher-order mode. -## Specialisations +The `computeHessian()` method takes a set of variables packaged in a `std::vector` and a function with signature `T foo(std::vector)`. -Hessians can be computed in `fwd_adj` or `fwd_fwd` higher-order mode. +## Return Types -The `computeHessian()` method takes a set of variables packaged in a `std::vector` and a function in the format `T foo(std::vector)`. -If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. +If provided with `RowIterators`, `computeHessian()` will write directly to them and return `void`. If no `RowIterators` are provided, the Hessian will be written to a `std::vector>` and returned. + +## Specialisations #### `fwd_adj` +```c++ +template +void computeHessian( + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + RowIterator first, RowIterator last, + xad::Tape> *tape = xad::Tape>::getActive()) +``` This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` -This mode does not require a Tape and can help reduce the overhead that comes with one. +```c++ +template +void computeHessian( + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + RowIterator first, RowIterator last) +``` +This mode does not require a Tape which can help reduce the overhead that comes with one. + + +## Example Use +Given $f(x, y, z, w) = sin(x * y) - cos(y * z) - sin(z * w) - cos(w * x)$, or +```c++ +auto foo = [](std::vector &x) -> AD +{ + return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); +}; +``` +with the derivatives of interest being +```c++ +std::vector x_ad({1.0, 1.5, 1.3, 1.2}); +``` +we'd like to compute the Hessian +$$ +H = \begin{bmatrix} +\frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial x \partial w} \\ +\frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^2 f}{\partial y^2} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial y \partial w} \\ +\frac{\partial^2 f}{\partial z \partial x} & \frac{\partial^2 f}{\partial z \partial y} & \frac{\partial^2 f}{\partial z^2} & \frac{\partial^2 f}{\partial z \partial w} \\ +\frac{\partial^2 f}{\partial w \partial x} & \frac{\partial^2 f}{\partial w \partial y} & \frac{\partial^2 f}{\partial w \partial z} & \frac{\partial^2 f}{\partial w^2} +\end{bmatrix} +$$ +First step is to setup the tape and active data types +```c++ + typedef xad::fwd_adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; +``` +Note that if no tape is setup, one will be created when computing the Hessian. `fwd_fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeHessian()`: +```c++ + auto foo = [](std::vector &x) -> AD + { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; + + auto hessian = xad::computeHessian(x_ad, foo); +``` +Note the signature of `foo()`. Any other signature will throw an error. + +This computes the relevant matrix +$$ +H = \begin{bmatrix} +-1.72257 & -1.42551 & 0 & 1.36687 \\ +-1.42551 & -1.6231 & 0.207107 & 0 \\ +0 & 0.207107 & 0.607009 & 1.54911 \\ +1.36687 & 0 & 1.54911 & 2.05226 +\end{bmatrix} +$$ +and prints it +```c++ + for (auto row : hessian) + { + for (auto elem : row) std::cout << elem << " "; + std::cout << std::endl; + } +``` diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index d443641..409014a 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -2,23 +2,89 @@ ## Overview -XAD implements a set of methods to compute the Jacobian matrix of a function and its inputs in `XAD/Jacobian.hpp`. +XAD implements a set of methods to compute the Jacobian matrix of a function in `XAD/Jacobian.hpp`. Note that the Jacobian header is not automatically included with `XAD/XAD.hpp`. Users must include it as needed. +Jacobians can be computed in `adj` or `fwd` mode. -## Specialisations +The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function with signature `std::vector foo(std::vector)`. -Jacobians can be computed in `adj` or `fwd` mode. +## Return Types -The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function in the format `std::vector foo(std::vector)`. -If provided with `RowIterator`s, `computeHessian()` will write directly to them and return `void`. If no `RowIterator`s are provided, the Hessian will be written to a `std::vector>` and returned. +If provided with `RowIterators`, `computeHessian()` will write directly to them and return `void`. If no `RowIterators` are provided, the Hessian will be written to a `std::vector>` and returned. + +## Specialisations -#### `adj` +#### `adj` +```c++ +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last, + xad::Tape *tape = xad::Tape::getActive()) +``` This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` +```c++ +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last) +``` This mode does not require a Tape and can help reduce the overhead that comes with one. + +## Example Use +Given $f(x, y, z, w) = [sin(x + y) sin(y + z) cos(z + w) cos(w + x)]$, or +```c++ +auto foo = [](std::vector &x) -> std::vector +{ + return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; +}; +``` +with the derivatives of interest being +```c++ +std::vector x_ad({1.0, 1.5, 1.3, 1.2}); +``` +we'd like to compute the Jacobian +$$ +J = \begin{bmatrix} \frac{\partial \sin(x + y)}{\partial x} & \frac{\partial \sin(x + y)}{\partial y} & \frac{\partial \sin(x + y)}{\partial z} & \frac{\partial \sin(x + y)}{\partial w} \\ \frac{\partial \sin(y + z)}{\partial x} & \frac{\partial \sin(y + z)}{\partial y} & \frac{\partial \sin(y + z)}{\partial z} & \frac{\partial \sin(y + z)}{\partial w} \\ \frac{\partial \cos(z + w)}{\partial x} & \frac{\partial \cos(z + w)}{\partial y} & \frac{\partial \cos(z + w)}{\partial z} & \frac{\partial \cos(z + w)}{\partial w} \\ \frac{\partial \cos(w + x)}{\partial x} & \frac{\partial \cos(w + x)}{\partial y} & \frac{\partial \cos(w + x)}{\partial z} & \frac{\partial \cos(w + x)}{\partial w} \end{bmatrix} +$$ +First step is to setup the tape and active data types +```c++ + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + tape_type tape; +``` +Note that if no tape is setup, one will be created when computing the Jacobian. `fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeJacobian()`: +```c++ + auto foo = [](std::vector &x) -> std::vector + { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; + + auto jacobian = xad::computeJacobian(x_ad, foo); +``` +Note the signature of `foo()`. Any other signature will throw an error. + +This computes the relevant matrix +$$ +\begin{bmatrix} +1 & 0.0707372 & 0 & 0 \\ +0 & 1 & 0.267499 & 0 \\ +0 & 0 & 1 & 0.362358 \\ +0.540302 & 0 & 0 & 1 +\end{bmatrix} +$$ +and prints it +```c++ + for (auto row : jacobian) + { + for (auto elem : row) std::cout << elem << " "; + std::cout << std::endl; + } +``` diff --git a/samples/CMakeLists.txt b/samples/CMakeLists.txt index 4511664..25076a2 100644 --- a/samples/CMakeLists.txt +++ b/samples/CMakeLists.txt @@ -70,5 +70,7 @@ add_subdirectory(checkpointing) add_subdirectory(external_function) add_subdirectory(SwapPricer) add_subdirectory(fwd_adj_2nd) +add_subdirectory(Hessian) +add_subdirectory(Jacobian) diff --git a/samples/Hessian/CMakeLists.txt b/samples/Hessian/CMakeLists.txt new file mode 100644 index 0000000..a4fb458 --- /dev/null +++ b/samples/Hessian/CMakeLists.txt @@ -0,0 +1,25 @@ +############################################################################## +# +# Hessian sample CMake file +# +# This file is part of XAD, a comprehensive C++ library for +# automatic differentiation. +# +# Copyright (C) 2010-2024 Xcelerit Computing Ltd. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero 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 Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +# +############################################################################## + +xad_add_sample(Hessian main.cpp) diff --git a/samples/Hessian/main.cpp b/samples/Hessian/main.cpp new file mode 100644 index 0000000..c7e1d1e --- /dev/null +++ b/samples/Hessian/main.cpp @@ -0,0 +1,59 @@ +/******************************************************************************* + + Computes + f(x, y, z, w) = sin(x * y) - cos(y * z) - sin(z * w) - cos(w * x) + and its Hessian matrix (second order derivatives) + using forward adjoint mode. + + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +#include +#include + +int main() +{ + // tape and active data type for 2nd order forward adjoint computation + typedef xad::fwd_adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + // initialize tape (not required) + tape_type tape; + + // input vector + std::vector x_ad({1.0, 1.5, 1.3, 1.2}); + + // many-input function + auto foo = [](std::vector &x) -> AD + { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; + + auto hessian = xad::computeHessian(x_ad, foo); + + // output results + for (auto row : hessian) + { + for (auto elem : row) std::cout << elem << " "; + std::cout << std::endl; + } +} diff --git a/samples/Jacobian/CMakeLists.txt b/samples/Jacobian/CMakeLists.txt new file mode 100644 index 0000000..471d5a4 --- /dev/null +++ b/samples/Jacobian/CMakeLists.txt @@ -0,0 +1,25 @@ +############################################################################## +# +# Jacobian sample CMake file +# +# This file is part of XAD, a comprehensive C++ library for +# automatic differentiation. +# +# Copyright (C) 2010-2024 Xcelerit Computing Ltd. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU Affero 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 Affero General Public License for more details. +# +# You should have received a copy of the GNU Affero General Public License +# along with this program. If not, see . +# +############################################################################## + +xad_add_sample(Jacobian main.cpp) diff --git a/samples/Jacobian/main.cpp b/samples/Jacobian/main.cpp new file mode 100644 index 0000000..df57967 --- /dev/null +++ b/samples/Jacobian/main.cpp @@ -0,0 +1,58 @@ +/******************************************************************************* + + Computes + f(x, y, z, w) = [sin(x + y) sin(y + z) cos(z + w) cos(w + x)] + and its Jacobian matrix using adjoint mode. + + + This file is part of XAD, a comprehensive C++ library for + automatic differentiation. + + Copyright (C) 2010-2024 Xcelerit Computing Ltd. + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU Affero 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 Affero General Public License for more details. + + You should have received a copy of the GNU Affero General Public License + along with this program. If not, see . + +******************************************************************************/ + +#include + +#include +#include + +int main() +{ + // tape and active data type for 1st order adjoint computation + typedef xad::adj mode; + typedef mode::tape_type tape_type; + typedef mode::active_type AD; + + // initialize tape (not required) + tape_type tape; + + // input vector + std::vector x_ad({1.0, 1.5, 1.3, 1.2}); + + // many-input many-output function + auto foo = [](std::vector &x) -> std::vector + { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; + + auto jacobian = xad::computeJacobian(x_ad, foo); + + // output results + for (auto row : jacobian) + { + for (auto elem : row) std::cout << elem << " "; + std::cout << std::endl; + } +} diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 0667577..4ae5972 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -58,6 +58,8 @@ void computeHessian( if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); + else if (std::distance(first->cbegin(), first->cend()) != domain) + throw OutOfRange("Iterator not allocated enough space"); static_assert( xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); @@ -112,6 +114,8 @@ void computeHessian( if (std::distance(first, last) != domain) throw OutOfRange("Iterator not allocated enough space"); + else if (std::distance(first->cbegin(), first->cend()) != domain) + throw OutOfRange("Iterator not allocated enough space"); static_assert( xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index bcaca9c..d26ff4e 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -53,12 +53,8 @@ void computeJacobian(const std::vector> &vec, RowIterator first, RowIterator last, xad::Tape *tape = xad::Tape::getActive()) { - auto v(vec); - unsigned int domain = static_cast(vec.size()), - codomain = static_cast(foo(v).size()); - - if (std::distance(first, last) != codomain) - throw OutOfRange("Iterator not allocated enough space"); + if (std::distance(first->cbegin(), first->cend()) != vec.size()) + throw OutOfRange("Iterator not allocated enough space (domain)"); static_assert( xad::detail::has_begin::value_type>::value, "RowIterator must dereference to a type that implements a begin() method"); @@ -69,21 +65,26 @@ void computeJacobian(const std::vector> &vec, tape = t.get(); } + auto v(vec); + unsigned int domain = static_cast(vec.size()), + codomain = static_cast(foo(v).size()); + + if (std::distance(first, last) != codomain) + throw OutOfRange("Iterator not allocated enough space (codomain)"); + tape->registerInputs(v); + tape->newRecording(); + auto y = foo(v); + tape->registerOutputs(y); auto row = first; for (unsigned int i = 0; i < codomain; i++, row++) { auto col = row->begin(); - for (unsigned int j = 0; j < domain; j++, col++) - { - tape->newRecording(); - xad::AReal y = foo(v)[i]; - tape->registerOutput(y); - derivative(y) = 1.0; - tape->computeAdjoints(); - *col = derivative(v[j]); - } + derivative(y[i]) = 1.0; + tape->computeAdjoints(); + for (unsigned int j = 0; j < domain; j++, col++) *col = derivative(v[j]); + tape->clearDerivatives(); } } @@ -105,15 +106,18 @@ void computeJacobian(const std::vector> &vec, std::function>(std::vector> &)> foo, RowIterator first, RowIterator last) { + if (std::distance(first->cbegin(), first->cend()) != vec.size()) + throw OutOfRange("Iterator not allocated enough space (domain)"); + static_assert( + xad::detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + auto v(vec); unsigned int domain = static_cast(vec.size()), codomain = static_cast(foo(v).size()); if (std::distance(first, last) != codomain) - throw OutOfRange("Iterator not allocated enough space"); - static_assert( - xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + throw OutOfRange("Iterator not allocated enough space (codomain)"); auto row = first; for (unsigned int i = 0; i < codomain; i++, row++) @@ -122,8 +126,9 @@ void computeJacobian(const std::vector> &vec, for (unsigned int j = 0; j < domain; j++, col++) { derivative(v[j]) = 1.0; - xad::FReal y = foo(v)[i]; - *col = derivative(y); + auto y = foo(v); + *col = derivative(y[i]); + std::cout << derivative(y[i]) << std::endl; derivative(v[j]) = 0.0; } } diff --git a/test/TypeTraits_test.cpp b/test/TypeTraits_test.cpp index c45af2d..617a268 100644 --- a/test/TypeTraits_test.cpp +++ b/test/TypeTraits_test.cpp @@ -28,10 +28,11 @@ #include #include -TEST(TypeTraits, HasBeginVectorOfVectorsIterator) +TEST(TypeTraits, HasBeginVectorOfVectorsIteratorDereferenced) { - using It = std::vector>::iterator; - EXPECT_FALSE(xad::detail::has_begin::value); + std::vector> t(0); + auto row = t.begin(); + EXPECT_TRUE(xad::detail::has_begin::value); } TEST(TypeTraits, HasBeginVectorOfVectors) @@ -40,7 +41,7 @@ TEST(TypeTraits, HasBeginVectorOfVectors) EXPECT_TRUE(xad::detail::has_begin::value); } -TEST(TypeTraits, HasBeginListIterator) +TEST(TypeTraits, HasBeginListIteratorReferenced) { using It = std::list::iterator; EXPECT_FALSE(xad::detail::has_begin::value); From 3d8be2ba69f4570b48a8e80d6bbdb40b49ebe14a Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 14:28:25 +0100 Subject: [PATCH 21/30] Added TypeTraits.hpp to CMakeLists --- src/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7b5321e..51d3fe4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -52,6 +52,7 @@ set(public_headers XAD/Tape.hpp XAD/TapeContainer.hpp XAD/Traits.hpp + XAD/TypeTraits.hpp XAD/UnaryExpr.hpp XAD/UnaryFunctors.hpp XAD/UnaryMathFunctors.hpp From 60f4245f52d09e1665151e5b555a307f8d19bebc Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 15:33:13 +0100 Subject: [PATCH 22/30] New Jacobian test, fixed Jacobian being made transposed and format changes to doc --- docs/ref/hessian.md | 16 +++++++++++++ docs/ref/jacobian.md | 18 ++++++++++++++- src/XAD/Jacobian.hpp | 19 ++++++++------- test/Jacobian_test.cpp | 52 ++++++++++++++++++++++++++++++++++-------- 4 files changed, 86 insertions(+), 19 deletions(-) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 89761d8..28694d4 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -27,6 +27,7 @@ void computeHessian( RowIterator first, RowIterator last, xad::Tape> *tape = xad::Tape>::getActive()) ``` + This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` @@ -38,22 +39,29 @@ void computeHessian( std::function>(std::vector>> &)> foo, RowIterator first, RowIterator last) ``` + This mode does not require a Tape which can help reduce the overhead that comes with one. ## Example Use + Given $f(x, y, z, w) = sin(x * y) - cos(y * z) - sin(z * w) - cos(w * x)$, or + ```c++ auto foo = [](std::vector &x) -> AD { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; ``` + with the derivatives of interest being + ```c++ std::vector x_ad({1.0, 1.5, 1.3, 1.2}); ``` + we'd like to compute the Hessian + $$ H = \begin{bmatrix} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial x \partial w} \\ @@ -62,7 +70,9 @@ H = \begin{bmatrix} \frac{\partial^2 f}{\partial w \partial x} & \frac{\partial^2 f}{\partial w \partial y} & \frac{\partial^2 f}{\partial w \partial z} & \frac{\partial^2 f}{\partial w^2} \end{bmatrix} $$ + First step is to setup the tape and active data types + ```c++ typedef xad::fwd_adj mode; typedef mode::tape_type tape_type; @@ -70,16 +80,20 @@ First step is to setup the tape and active data types tape_type tape; ``` + Note that if no tape is setup, one will be created when computing the Hessian. `fwd_fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeHessian()`: + ```c++ auto foo = [](std::vector &x) -> AD { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; auto hessian = xad::computeHessian(x_ad, foo); ``` + Note the signature of `foo()`. Any other signature will throw an error. This computes the relevant matrix + $$ H = \begin{bmatrix} -1.72257 & -1.42551 & 0 & 1.36687 \\ @@ -88,7 +102,9 @@ H = \begin{bmatrix} 1.36687 & 0 & 1.54911 & 2.05226 \end{bmatrix} $$ + and prints it + ```c++ for (auto row : hessian) { diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 409014a..3583f88 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -17,8 +17,8 @@ If provided with `RowIterators`, `computeHessian()` will write directly to them ## Specialisations - #### `adj` + ```c++ template void computeJacobian(const std::vector> &vec, @@ -26,6 +26,7 @@ void computeJacobian(const std::vector> &vec, RowIterator first, RowIterator last, xad::Tape *tape = xad::Tape::getActive()) ``` + This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. #### `fwd_fwd` @@ -36,25 +37,34 @@ void computeJacobian(const std::vector> &vec, std::function>(std::vector> &)> foo, RowIterator first, RowIterator last) ``` + This mode does not require a Tape and can help reduce the overhead that comes with one. ## Example Use + Given $f(x, y, z, w) = [sin(x + y) sin(y + z) cos(z + w) cos(w + x)]$, or + ```c++ auto foo = [](std::vector &x) -> std::vector { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; ``` + with the derivatives of interest being + ```c++ std::vector x_ad({1.0, 1.5, 1.3, 1.2}); ``` + we'd like to compute the Jacobian + $$ J = \begin{bmatrix} \frac{\partial \sin(x + y)}{\partial x} & \frac{\partial \sin(x + y)}{\partial y} & \frac{\partial \sin(x + y)}{\partial z} & \frac{\partial \sin(x + y)}{\partial w} \\ \frac{\partial \sin(y + z)}{\partial x} & \frac{\partial \sin(y + z)}{\partial y} & \frac{\partial \sin(y + z)}{\partial z} & \frac{\partial \sin(y + z)}{\partial w} \\ \frac{\partial \cos(z + w)}{\partial x} & \frac{\partial \cos(z + w)}{\partial y} & \frac{\partial \cos(z + w)}{\partial z} & \frac{\partial \cos(z + w)}{\partial w} \\ \frac{\partial \cos(w + x)}{\partial x} & \frac{\partial \cos(w + x)}{\partial y} & \frac{\partial \cos(w + x)}{\partial z} & \frac{\partial \cos(w + x)}{\partial w} \end{bmatrix} $$ + First step is to setup the tape and active data types + ```c++ typedef xad::adj mode; typedef mode::tape_type tape_type; @@ -62,16 +72,20 @@ First step is to setup the tape and active data types tape_type tape; ``` + Note that if no tape is setup, one will be created when computing the Jacobian. `fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeJacobian()`: + ```c++ auto foo = [](std::vector &x) -> std::vector { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; auto jacobian = xad::computeJacobian(x_ad, foo); ``` + Note the signature of `foo()`. Any other signature will throw an error. This computes the relevant matrix + $$ \begin{bmatrix} 1 & 0.0707372 & 0 & 0 \\ @@ -80,7 +94,9 @@ $$ 0.540302 & 0 & 0 & 1 \end{bmatrix} $$ + and prints it + ```c++ for (auto row : jacobian) { diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index d26ff4e..5b4debb 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -28,6 +28,7 @@ #include #include +#include #include #include @@ -120,16 +121,18 @@ void computeJacobian(const std::vector> &vec, throw OutOfRange("Iterator not allocated enough space (codomain)"); auto row = first; - for (unsigned int i = 0; i < codomain; i++, row++) + for (unsigned int i = 0; i < domain; i++) { - auto col = row->begin(); - for (unsigned int j = 0; j < domain; j++, col++) + derivative(v[i]) = 1.0; + auto y = foo(v); + derivative(v[i]) = 0.0; + for (unsigned int j = 0; j < codomain; j++) { - derivative(v[j]) = 1.0; - auto y = foo(v); - *col = derivative(y[i]); - std::cout << derivative(y[i]) << std::endl; - derivative(v[j]) = 0.0; + row = first; + std::advance(row, j); + auto col = row->begin(); + std::advance(col, i); + *col = derivative(y[j]); } } } diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 323fb3f..b9dab9c 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -47,7 +47,8 @@ TEST(JacobianTest, SimpleAdjoint) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::vector> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + std::vector> expected_jacobian = {{1.0, cos(x[1])}, // + {cos(x[0]), 1.0}}; auto computed_jacobian = xad::computeJacobian(x, foo, &tape); @@ -70,7 +71,8 @@ TEST(JacobianTest, SimpleAdjointIteratorAutoTape) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::list> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + std::list> expected_jacobian = {{1.0, cos(x[1])}, // + {cos(x[0]), 1.0}}; std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( @@ -102,7 +104,8 @@ TEST(JacobianTest, SimpleForward) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::vector> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + std::vector> expected_jacobian = {{1.0, cos(x[1])}, // + {cos(x[0]), 1.0}}; auto computed_jacobian = xad::computeJacobian(x, foo); @@ -122,7 +125,8 @@ TEST(JacobianTest, SimpleForwardIterator) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::list> expected_jacobian = {{1.0, cos(x[1])}, {cos(x[0]), 1.0}}; + std::list> expected_jacobian = {{1.0, cos(x[1])}, // + {cos(x[0]), 1.0}}; std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( @@ -177,7 +181,28 @@ TEST(JacobianTest, DomainLargerThanCodomainForward) // f(x) = [ x[0] + x[1], x[2] * x[3] ] auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[2] * x[3]}; }; - std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, {0.0, 0.0, x[3], x[2]}}; + std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // + {0.0, 0.0, x[3], x[2]}}; + + auto computed_jacobian = xad::computeJacobian(x, foo); + + for (unsigned int i = 0; i < expected_jacobian.size(); i++) + for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); +} + +TEST(JacobianTest, DomainLargerThanCodomainAdjoint) +{ + typedef xad::adj mode; + typedef mode::active_type AD; + + std::vector x = {1.0, 2.0, 3.0, 4.0}; + + // f(x) = [ x[0] + x[1], x[2] * x[3] ] + auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[2] * x[3]}; }; + + std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // + {0.0, 0.0, x[3], x[2]}}; auto computed_jacobian = xad::computeJacobian(x, foo); @@ -200,7 +225,9 @@ TEST(JacobianTest, DomainSmallerThanCodomainAdjoint) auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[0] - x[1], x[0] * x[1]}; }; - std::vector> expected_jacobian = {{1.0, 1.0}, {1.0, -1.0}, {x[1], x[0]}}; + std::vector> expected_jacobian = {{1.0, 1.0}, // + {1.0, -1.0}, + {x[1], x[0]}}; auto computed_jacobian = xad::computeJacobian(x, foo, &tape); @@ -220,8 +247,9 @@ TEST(JacobianTest, ComplexDomainNotEqualCodomainForwardIterator) auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[1] * x[2], exp(x[0])}; }; - std::list> expected_jacobian = { - {1.0, 1.0, 0.0}, {0.0, x[2], x[1]}, {exp(x[0]), 0.0, 0.0}}; + std::list> expected_jacobian = {{1.0, 1.0, 0.0}, // + {0.0, x[2], x[1]}, + {exp(x[0]), 0.0, 0.0}}; std::list> computed_jacobian(3, std::list(3, 0.0)); xad::computeJacobian( @@ -252,7 +280,8 @@ TEST(JacobianTest, TrigonometricFunctionForward) // f(x) = [ sin(x[0]), cos(x[1]) ] auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; - std::vector> expected_jacobian = {{cos(x[0]), 0.0}, {0.0, -sin(x[1])}}; + std::vector> expected_jacobian = {{cos(x[0]), 0.0}, // + {0.0, -sin(x[1])}}; auto computed_jacobian = xad::computeJacobian(x, foo); @@ -271,11 +300,14 @@ TEST(JacobianTest, TrigonometricFunctionAdjointAutoTape) // f(x) = [ sin(x[0]), cos(x[1]) ] auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; - std::vector> expected_jacobian = {{cos(x[0]), 0.0}, {0.0, -sin(x[1])}}; + std::vector> expected_jacobian = {{cos(x[0]), 0.0}, // + {0.0, -sin(x[1])}}; auto computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) + { ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + } } From e5875bdcd897a2ef965c90f7d95eece232ddf57a Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 15:36:05 +0100 Subject: [PATCH 23/30] Fixed display of Jacobian in doc --- docs/ref/jacobian.md | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 3583f88..2d01592 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -60,7 +60,12 @@ std::vector x_ad({1.0, 1.5, 1.3, 1.2}); we'd like to compute the Jacobian $$ -J = \begin{bmatrix} \frac{\partial \sin(x + y)}{\partial x} & \frac{\partial \sin(x + y)}{\partial y} & \frac{\partial \sin(x + y)}{\partial z} & \frac{\partial \sin(x + y)}{\partial w} \\ \frac{\partial \sin(y + z)}{\partial x} & \frac{\partial \sin(y + z)}{\partial y} & \frac{\partial \sin(y + z)}{\partial z} & \frac{\partial \sin(y + z)}{\partial w} \\ \frac{\partial \cos(z + w)}{\partial x} & \frac{\partial \cos(z + w)}{\partial y} & \frac{\partial \cos(z + w)}{\partial z} & \frac{\partial \cos(z + w)}{\partial w} \\ \frac{\partial \cos(w + x)}{\partial x} & \frac{\partial \cos(w + x)}{\partial y} & \frac{\partial \cos(w + x)}{\partial z} & \frac{\partial \cos(w + x)}{\partial w} \end{bmatrix} +J = \begin{bmatrix} +\frac{\partial \sin(x + y)}{\partial x} & \frac{\partial \sin(x + y)}{\partial y} & \frac{\partial \sin(x + y)}{\partial z} & \frac{\partial \sin(x + y)}{\partial w} \\ +\frac{\partial \sin(y + z)}{\partial x} & \frac{\partial \sin(y + z)}{\partial y} & \frac{\partial \sin(y + z)}{\partial z} & \frac{\partial \sin(y + z)}{\partial w} \\ +\frac{\partial \cos(z + w)}{\partial x} & \frac{\partial \cos(z + w)}{\partial y} & \frac{\partial \cos(z + w)}{\partial z} & \frac{\partial \cos(z + w)}{\partial w} \\ +\frac{\partial \cos(w + x)}{\partial x} & \frac{\partial \cos(w + x)}{\partial y} & \frac{\partial \cos(w + x)}{\partial z} & \frac{\partial \cos(w + x)}{\partial w} +\end{bmatrix} $$ First step is to setup the tape and active data types From fdf7447350f7abba57a5b7f3bb9f7a805631e401 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 16:10:14 +0100 Subject: [PATCH 24/30] Added tests for OutOfBounds exceptions --- test/Hessian_test.cpp | 25 ++++++++++++++++++++++ test/Jacobian_test.cpp | 48 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 12aa2d8..f67ee86 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -384,3 +384,28 @@ TEST(HessianTest, QuadraticForwardAdjointAutoTape) for (unsigned int j = 0; j < expected_hessian[i].size(); j++) ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } + +TEST(HessianTest, OutOfBoundsDomainSizeMismatch) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + + std::vector x = {1.0, 2.0}; + + auto foo = [](std::vector &x) -> AD { return x[0]; }; + + std::vector> jacobian(2, std::vector(3)); + + auto launch = [](std::vector x, + std::function>( + std::vector>> &)> + foo, + std::vector>>>::iterator first, + std::vector>>>::iterator last) + { + using RowIterator = decltype(first); + xad::computeHessian(x, foo, first, last); + }; + + EXPECT_THROW(launch(x, foo, begin(jacobian), end(jacobian)), xad::OutOfRange); +} diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index b9dab9c..f9cffd8 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -311,3 +311,51 @@ TEST(JacobianTest, TrigonometricFunctionAdjointAutoTape) ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); } } + +TEST(JacobianTest, OutOfBoundsDomainSizeMismatch) +{ + typedef xad::adj mode; + typedef mode::active_type AD; + + std::vector x = {1.0, 2.0}; + + auto foo = [](std::vector &x) -> std::vector { return {x[0], x[1]}; }; + + std::vector> jacobian(2, std::vector(3)); + + auto launch = + [](std::vector x, + std::function>(std::vector> &)> foo, + std::vector>>::iterator first, + std::vector>>::iterator last) + { + using RowIterator = decltype(first); + xad::computeJacobian(x, foo, first, last); + }; + + EXPECT_THROW(launch(x, foo, begin(jacobian), end(jacobian)), xad::OutOfRange); +} + +TEST(JacobianTest, OutOfBoundsCodomainSizeMismatch) +{ + typedef xad::adj mode; + typedef mode::active_type AD; + + std::vector x = {1.0}; + + auto foo = [](std::vector &x) -> std::vector { return {x[0], x[0]}; }; + + std::vector> jacobian(1, std::vector(1)); + + auto launch = + [](std::vector x, + std::function>(std::vector> &)> foo, + std::vector>>::iterator first, + std::vector>>::iterator last) + { + using RowIterator = decltype(first); + xad::computeJacobian(x, foo, first, last); + }; + + EXPECT_THROW(launch(x, foo, begin(jacobian), end(jacobian)), xad::OutOfRange); +} From 7fa8fa1e19d71cf9a138d1293b9ec9092d8b0baf Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 16:25:46 +0100 Subject: [PATCH 25/30] More exhaustive tests --- src/XAD/Hessian.hpp | 10 +++---- test/Hessian_test.cpp | 60 +++++++++++++++++++++++++++++------------- test/Jacobian_test.cpp | 35 +++++++++++++++++++++++- 3 files changed, 79 insertions(+), 26 deletions(-) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 4ae5972..4888622 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -56,9 +56,8 @@ void computeHessian( { unsigned int domain(static_cast(vec.size())); - if (std::distance(first, last) != domain) - throw OutOfRange("Iterator not allocated enough space"); - else if (std::distance(first->cbegin(), first->cend()) != domain) + if (std::distance(first, last) != domain || + std::distance(first->cbegin(), first->cend()) != domain) throw OutOfRange("Iterator not allocated enough space"); static_assert( xad::detail::has_begin::value_type>::value, @@ -112,9 +111,8 @@ void computeHessian( { unsigned int domain(static_cast(vec.size())); - if (std::distance(first, last) != domain) - throw OutOfRange("Iterator not allocated enough space"); - else if (std::distance(first->cbegin(), first->cend()) != domain) + if (std::distance(first, last) != domain || + std::distance(first->cbegin(), first->cend()) != domain) throw OutOfRange("Iterator not allocated enough space"); static_assert( xad::detail::has_begin::value_type>::value, diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index f67ee86..0749db7 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -58,6 +58,47 @@ TEST(HessianTest, QuadraticForwardAdjoint) ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } +TEST(HessianTest, QuadraticForwardAdjointAutoTape) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; + + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + auto computed_hessian = xad::computeHessian(x, foo); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + +TEST(HessianTest, QuadraticForwardAdjointFetchTape) +{ + typedef xad::fwd_adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + + tape_type tape; + + std::vector x = {3.0, 2.0}; + + // f(x) = x[0]^2 + x[1]^2 + auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; + + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + + auto computed_hessian = xad::computeHessian(x, foo); + + for (unsigned int i = 0; i < expected_hessian.size(); i++) + for (unsigned int j = 0; j < expected_hessian[i].size(); j++) + ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); +} + TEST(HessianTest, QuadraticForwardAdjointWithIterator) { typedef xad::fwd_adj mode; @@ -366,25 +407,6 @@ TEST(HessianTest, LargeHessianForwardForward) ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); } -TEST(HessianTest, QuadraticForwardAdjointAutoTape) -{ - typedef xad::fwd_adj mode; - typedef mode::active_type AD; - - std::vector x = {3.0, 2.0}; - - // f(x) = x[0]^2 + x[1]^2 - auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - - std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - - auto computed_hessian = xad::computeHessian(x, foo); - - for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); -} - TEST(HessianTest, OutOfBoundsDomainSizeMismatch) { typedef xad::fwd_adj mode; diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index f9cffd8..12bcca2 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -60,9 +60,42 @@ TEST(JacobianTest, SimpleAdjoint) TEST(JacobianTest, SimpleAdjointIteratorAutoTape) { typedef xad::adj mode; - typedef mode::tape_type tape_type; typedef mode::active_type AD; + std::vector x = {3.0, 1.0}; + + // f(x) = [ x[0] + sin(x[1]), x[1] + sin(x[0]) ] + auto foo = [](std::vector &x) -> std::vector + { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; + + std::list> expected_jacobian = {{1.0, cos(x[1])}, // + {cos(x[0]), 1.0}}; + + std::list> computed_jacobian(2, std::list(2, 0.0)); + xad::computeJacobian( + x, foo, begin(computed_jacobian), end(computed_jacobian)); + + auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); + while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) + { + auto col1 = row1->begin(), col2 = row2->begin(); + while (col1 != row1->end() && col2 != row2->end()) + { + ASSERT_EQ(*col1, *col2); + col1++; + col2++; + } + row1++; + row2++; + } +} + +TEST(JacobianTest, SimpleAdjointIteratorFetchTape) +{ + typedef xad::adj mode; + typedef mode::active_type AD; + typedef mode::tape_type tape_type; + tape_type tape; std::vector x = {3.0, 1.0}; From 7c0382ce27f0a2a95e8e93f3699c966202849d92 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 21:41:56 +0100 Subject: [PATCH 26/30] Minimised total function calls in Jacobian at expense of late iterator distance check --- src/XAD/Jacobian.hpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 5b4debb..1e66f3f 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -67,15 +67,14 @@ void computeJacobian(const std::vector> &vec, } auto v(vec); - unsigned int domain = static_cast(vec.size()), - codomain = static_cast(foo(v).size()); - - if (std::distance(first, last) != codomain) - throw OutOfRange("Iterator not allocated enough space (codomain)"); tape->registerInputs(v); tape->newRecording(); auto y = foo(v); + unsigned int domain = static_cast(vec.size()), + codomain = static_cast(y.size()); + if (std::distance(first, last) != codomain) + throw OutOfRange("Iterator not allocated enough space (codomain)"); tape->registerOutputs(y); auto row = first; From 06ee32fae97da0310ccdb18eebd8db2b527543c5 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 22:08:11 +0100 Subject: [PATCH 27/30] Re-format as per Codacy --- docs/ref/headers.md | 23 ++++++++++---------- docs/ref/hessian.md | 50 +++++++++++++++++++++++++++++++------------- docs/ref/jacobian.md | 48 +++++++++++++++++++++++++++++++----------- 3 files changed, 83 insertions(+), 38 deletions(-) diff --git a/docs/ref/headers.md b/docs/ref/headers.md index 112968b..4e51c09 100644 --- a/docs/ref/headers.md +++ b/docs/ref/headers.md @@ -9,14 +9,15 @@ Typically, this is all that clients need to include. There are four additional headers provided that can be included on demand: -* `XAD/Complex.hpp` - For using complex numbers with XAD data types (see [Complex](complex.md)). - This header should be included wherever [`#!c++ std::complex`](https://en.cppreference.com/w/cpp/numeric/complex) is used. -* `XAD/StdCompatibility.hpp` - This header imports the XAD math functions - into the `std` namespace, for compatibility reasons. - It enables using constructs like [`#!c++ std::sin(x)`](https://en.cppreference.com/w/cpp/numeric/math/sin) where `x` is an XAD type. - Additionally, it also specialises [`#!c++ std::numeric_limits`](https://en.cppreference.com/w/cpp/types/numeric_limits) for the XAD data types, - so that it provides traits similar to the standard floating point types. -* `XAD/Hessian.hpp` - Imports methods for computing the Hessian matrix of a single output function into the `xad` - namespace. -* `XAD/Jacobian.hpp` - Imports methods for computing the Jacobian matrix of a function with multiple inputs and multiple - outputs into the `xad` namespace. +* `XAD/Complex.hpp` - For using complex numbers with XAD data types + (see [Complex](complex.md)). + This header should be included wherever [`#!c++ std::complex`](https://en.cppreference.com/w/cpp/numeric/complex) is used. +* `XAD/StdCompatibility.hpp` - This header imports the XAD math functions + into the `std` namespace, for compatibility reasons. + It enables using constructs like [`#!c++ std::sin(x)`](https://en.cppreference.com/w/cpp/numeric/math/sin) where `x` is an XAD type. + Additionally, it also specialises [`#!c++ std::numeric_limits`](https://en.cppreference.com/w/cpp/types/numeric_limits) for the XAD data types, + so that it provides traits similar to the standard floating point types. +* `XAD/Hessian.hpp` - Imports methods for computing the Hessian matrix of a + single output function into the `xad` namespace. +* `XAD/Jacobian.hpp` - Imports methods for computing the Jacobian matrix of a + function with multiple inputs and multiple outputs into the `xad` namespace. diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 28694d4..c76d379 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -9,11 +9,14 @@ Users must include it as needed. Hessians can be computed in `fwd_adj` or `fwd_fwd` higher-order mode. -The `computeHessian()` method takes a set of variables packaged in a `std::vector` and a function with signature `T foo(std::vector)`. +The `computeHessian()` method takes a set of variables packaged in a +`std::vector` and a function with signature `T foo(std::vector)`. ## Return Types -If provided with `RowIterators`, `computeHessian()` will write directly to them and return `void`. If no `RowIterators` are provided, the Hessian will be written to a `std::vector>` and returned. +If provided with `RowIterators`, `computeHessian()` will write directly to +them and return `void`. If no `RowIterators` are provided, the Hessian will +be written to a `std::vector>` and returned. ## Specialisations @@ -22,26 +25,28 @@ If provided with `RowIterators`, `computeHessian()` will write directly to them ```c++ template void computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo, + const std::vector &vec, + std::function(std::vector &)> foo, RowIterator first, RowIterator last, xad::Tape> *tape = xad::Tape>::getActive()) ``` -This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. +This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape +will be instantiated within the method or set to the current active Tape using +`Tape::getActive()` if none is passed as argument. #### `fwd_fwd` ```c++ template void computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo, + const std::vector &vec, + std::function(std::vector &)> foo, RowIterator first, RowIterator last) ``` -This mode does not require a Tape which can help reduce the overhead that comes with one. - +This mode does not require a Tape which can help reduce the overhead that comes +with one. ## Example Use @@ -50,7 +55,8 @@ Given $f(x, y, z, w) = sin(x * y) - cos(y * z) - sin(z * w) - cos(w * x)$, or ```c++ auto foo = [](std::vector &x) -> AD { - return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); + return sin(x[0] * x[1]) - cos(x[1] * x[2]) + - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; ``` @@ -64,10 +70,22 @@ we'd like to compute the Hessian $$ H = \begin{bmatrix} -\frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} & \frac{\partial^2 f}{\partial x \partial z} & \frac{\partial^2 f}{\partial x \partial w} \\ -\frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^2 f}{\partial y^2} & \frac{\partial^2 f}{\partial y \partial z} & \frac{\partial^2 f}{\partial y \partial w} \\ -\frac{\partial^2 f}{\partial z \partial x} & \frac{\partial^2 f}{\partial z \partial y} & \frac{\partial^2 f}{\partial z^2} & \frac{\partial^2 f}{\partial z \partial w} \\ -\frac{\partial^2 f}{\partial w \partial x} & \frac{\partial^2 f}{\partial w \partial y} & \frac{\partial^2 f}{\partial w \partial z} & \frac{\partial^2 f}{\partial w^2} +\frac{\partial^2 f}{\partial x^2} & +\frac{\partial^2 f}{\partial x \partial y} & +\frac{\partial^2 f}{\partial x \partial z} & +\frac{\partial^2 f}{\partial x \partial w} \\ +\frac{\partial^2 f}{\partial y \partial x} & +\frac{\partial^2 f}{\partial y^2} & +\frac{\partial^2 f}{\partial y \partial z} & +\frac{\partial^2 f}{\partial y \partial w} \\ +\frac{\partial^2 f}{\partial z \partial x} & +\frac{\partial^2 f}{\partial z \partial y} & +\frac{\partial^2 f}{\partial z^2} & +\frac{\partial^2 f}{\partial z \partial w} \\ +\frac{\partial^2 f}{\partial w \partial x} & +\frac{\partial^2 f}{\partial w \partial y} & +\frac{\partial^2 f}{\partial w \partial z} & +\frac{\partial^2 f}{\partial w^2} \end{bmatrix} $$ @@ -81,7 +99,9 @@ First step is to setup the tape and active data types tape_type tape; ``` -Note that if no tape is setup, one will be created when computing the Hessian. `fwd_fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeHessian()`: +Note that if no tape is setup, one will be created when computing the Hessian. +`fwd_fwd` mode is also supported in the same fashion. All that is left to do +is define our input values and our function, then call `computeHessian()`: ```c++ auto foo = [](std::vector &x) -> AD diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 2d01592..cda367b 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -9,11 +9,15 @@ Users must include it as needed. Jacobians can be computed in `adj` or `fwd` mode. -The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function with signature `std::vector foo(std::vector)`. +The `computeJacobian()` method takes a set of variables packaged in a +`std::vector` and a function with signature +`std::vector foo(std::vector)`. ## Return Types -If provided with `RowIterators`, `computeHessian()` will write directly to them and return `void`. If no `RowIterators` are provided, the Hessian will be written to a `std::vector>` and returned. +If provided with `RowIterators`, `computeHessian()` will write directly to +them and return `void`. If no `RowIterators` are provided, the Hessian will be +written to a `std::vector>` and returned. ## Specialisations @@ -21,13 +25,15 @@ If provided with `RowIterators`, `computeHessian()` will write directly to them ```c++ template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, +void computeJacobian(const std::vector &vec, + std::function(std::vector &)> foo, RowIterator first, RowIterator last, xad::Tape *tape = xad::Tape::getActive()) ``` -This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. +This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will +be instantiated within the method or set to the current active Tape using +`Tape::getActive()` if none is passed as argument. #### `fwd_fwd` @@ -38,7 +44,8 @@ void computeJacobian(const std::vector> &vec, RowIterator first, RowIterator last) ``` -This mode does not require a Tape and can help reduce the overhead that comes with one. +This mode does not require a Tape and can help reduce the overhead that +comes with one. ## Example Use @@ -61,10 +68,22 @@ we'd like to compute the Jacobian $$ J = \begin{bmatrix} -\frac{\partial \sin(x + y)}{\partial x} & \frac{\partial \sin(x + y)}{\partial y} & \frac{\partial \sin(x + y)}{\partial z} & \frac{\partial \sin(x + y)}{\partial w} \\ -\frac{\partial \sin(y + z)}{\partial x} & \frac{\partial \sin(y + z)}{\partial y} & \frac{\partial \sin(y + z)}{\partial z} & \frac{\partial \sin(y + z)}{\partial w} \\ -\frac{\partial \cos(z + w)}{\partial x} & \frac{\partial \cos(z + w)}{\partial y} & \frac{\partial \cos(z + w)}{\partial z} & \frac{\partial \cos(z + w)}{\partial w} \\ -\frac{\partial \cos(w + x)}{\partial x} & \frac{\partial \cos(w + x)}{\partial y} & \frac{\partial \cos(w + x)}{\partial z} & \frac{\partial \cos(w + x)}{\partial w} +\frac{\partial \sin(x + y)}{\partial x} & +\frac{\partial \sin(x + y)}{\partial y} & +\frac{\partial \sin(x + y)}{\partial z} & +\frac{\partial \sin(x + y)}{\partial w} \\ +\frac{\partial \sin(y + z)}{\partial x} & +\frac{\partial \sin(y + z)}{\partial y} & +\frac{\partial \sin(y + z)}{\partial z} & +\frac{\partial \sin(y + z)}{\partial w} \\ +\frac{\partial \cos(z + w)}{\partial x} & +\frac{\partial \cos(z + w)}{\partial y} & +\frac{\partial \cos(z + w)}{\partial z} & +\frac{\partial \cos(z + w)}{\partial w} \\ +\frac{\partial \cos(w + x)}{\partial x} & +\frac{\partial \cos(w + x)}{\partial y} & +\frac{\partial \cos(w + x)}{\partial z} & +\frac{\partial \cos(w + x)}{\partial w} \end{bmatrix} $$ @@ -78,11 +97,16 @@ First step is to setup the tape and active data types tape_type tape; ``` -Note that if no tape is setup, one will be created when computing the Jacobian. `fwd` mode is also supported in the same fashion. All that is left to do is define our input values and our function, then call `computeJacobian()`: +Note that if no tape is setup, one will be created when computing the Jacobian. +`fwd` mode is also supported in the same fashion. All that is left to do is +define our input values and our function, then call `computeJacobian()`: ```c++ auto foo = [](std::vector &x) -> std::vector - { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; + { return {sin(x[0] + x[1]), + sin(x[1] + x[2]), + cos(x[2] + x[3]), + cos(x[3] + x[0])}; }; auto jacobian = xad::computeJacobian(x_ad, foo); ``` From 0cdeda2468e1e8036d9c46364c5c1d1b96d9dc1f Mon Sep 17 00:00:00 2001 From: imraneamri Date: Tue, 11 Jun 2024 22:22:36 +0100 Subject: [PATCH 28/30] More Markdown linting --- docs/README.md | 30 +++++++++++++++--------------- docs/ref/headers.md | 22 +++++++++++----------- docs/ref/hessian.md | 3 ++- docs/ref/jacobian.md | 4 ++-- 4 files changed, 30 insertions(+), 29 deletions(-) diff --git a/docs/README.md b/docs/README.md index 456ae70..e0fa74b 100644 --- a/docs/README.md +++ b/docs/README.md @@ -3,18 +3,18 @@ This folder contains the reference documentation for XAD, which is imported to the website at [auto-differentiation.github.io](https://auto-differentiation.github.io). -* [Headers and Namespaces](ref/headers.md) -* [AD Mode Interface](ref/interface.md) -* [Forward Mode Type FReal](ref/freal.md) -* [Adjoint Mode Type AReal](ref/areal.md) -* [Global Functions](ref/global.md) -* [Mathematical Operations](ref/math.md) -* [Complex](ref/complex.md) -* [Smoothed Mathematical Functions](ref/smooth-math.md) -* [Expressions](ref/expressions.md) -* [Tape](ref/tape.md) -* [CheckpointCallback](ref/chkpt_cb.md) -* [Exceptions](ref/exceptions.md) -* [Hessian](ref/hessian.md) -* [Jacobian](ref/jacobian.md) -* [Version Information](ref/version.md) +* [Headers and Namespaces](ref/headers.md) +* [AD Mode Interface](ref/interface.md) +* [Forward Mode Type FReal](ref/freal.md) +* [Adjoint Mode Type AReal](ref/areal.md) +* [Global Functions](ref/global.md) +* [Mathematical Operations](ref/math.md) +* [Complex](ref/complex.md) +* [Smoothed Mathematical Functions](ref/smooth-math.md) +* [Expressions](ref/expressions.md) +* [Tape](ref/tape.md) +* [CheckpointCallback](ref/chkpt_cb.md) +* [Exceptions](ref/exceptions.md) +* [Hessian](ref/hessian.md) +* [Jacobian](ref/jacobian.md) +* [Version Information](ref/version.md) diff --git a/docs/ref/headers.md b/docs/ref/headers.md index 4e51c09..acaf9c3 100644 --- a/docs/ref/headers.md +++ b/docs/ref/headers.md @@ -10,14 +10,14 @@ Typically, this is all that clients need to include. There are four additional headers provided that can be included on demand: * `XAD/Complex.hpp` - For using complex numbers with XAD data types - (see [Complex](complex.md)). - This header should be included wherever [`#!c++ std::complex`](https://en.cppreference.com/w/cpp/numeric/complex) is used. -* `XAD/StdCompatibility.hpp` - This header imports the XAD math functions - into the `std` namespace, for compatibility reasons. - It enables using constructs like [`#!c++ std::sin(x)`](https://en.cppreference.com/w/cpp/numeric/math/sin) where `x` is an XAD type. - Additionally, it also specialises [`#!c++ std::numeric_limits`](https://en.cppreference.com/w/cpp/types/numeric_limits) for the XAD data types, - so that it provides traits similar to the standard floating point types. -* `XAD/Hessian.hpp` - Imports methods for computing the Hessian matrix of a - single output function into the `xad` namespace. -* `XAD/Jacobian.hpp` - Imports methods for computing the Jacobian matrix of a - function with multiple inputs and multiple outputs into the `xad` namespace. + (see [Complex](complex.md)). + This header should be included wherever [`#!c++ std::complex`](https://en.cppreference.com/w/cpp/numeric/complex) is used. +* `XAD/StdCompatibility.hpp` - This header imports the XAD math functions + into the `std` namespace, for compatibility reasons. + It enables using constructs like [`#!c++ std::sin(x)`](https://en.cppreference.com/w/cpp/numeric/math/sin) where `x` is an XAD type. + Additionally, it also specialises [`#!c++ std::numeric_limits`](https://en.cppreference.com/w/cpp/types/numeric_limits) for the XAD data types, + so that it provides traits similar to the standard floating point types. +* `XAD/Hessian.hpp` - Imports methods for computing the Hessian matrix of a + single output function into the `xad` namespace. +* `XAD/Jacobian.hpp` - Imports methods for computing the Jacobian matrix of a + function with multiple inputs and multiple outputs into the `xad` namespace. diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index c76d379..4853266 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -105,7 +105,8 @@ is define our input values and our function, then call `computeHessian()`: ```c++ auto foo = [](std::vector &x) -> AD - { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; + { return sin(x[0] * x[1]) - cos(x[1] * x[2]) + - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; auto hessian = xad::computeHessian(x_ad, foo); ``` diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index cda367b..9006f9a 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -39,8 +39,8 @@ be instantiated within the method or set to the current active Tape using ```c++ template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, +void computeJacobian(const std::vector &vec, + std::function(std::vector &)> foo, RowIterator first, RowIterator last) ``` From 3fa8a7a9657a9e66f5546fe4afbbed08182a9897 Mon Sep 17 00:00:00 2001 From: imraneamri Date: Wed, 12 Jun 2024 10:40:27 +0100 Subject: [PATCH 29/30] Increased test rigor, Hessian & Jacobian final tweaks and doc fixes --- docs/ref/hessian.md | 32 ++++--- docs/ref/jacobian.md | 44 +++++---- samples/Hessian/main.cpp | 4 +- samples/Jacobian/main.cpp | 4 +- src/XAD/Hessian.hpp | 49 +++++----- src/XAD/Jacobian.hpp | 39 ++++---- test/Hessian_test.cpp | 168 ++++++++++++++------------------- test/Jacobian_test.cpp | 192 ++++++++++++++++---------------------- 8 files changed, 236 insertions(+), 296 deletions(-) diff --git a/docs/ref/hessian.md b/docs/ref/hessian.md index 4853266..6eff8bb 100644 --- a/docs/ref/hessian.md +++ b/docs/ref/hessian.md @@ -16,41 +16,43 @@ The `computeHessian()` method takes a set of variables packaged in a If provided with `RowIterators`, `computeHessian()` will write directly to them and return `void`. If no `RowIterators` are provided, the Hessian will -be written to a `std::vector>` and returned. +be written to a `std::vector>` and returned (`T` is +usually `double`). ## Specialisations -#### `fwd_adj` +### Forward over Adjoint Mode ```c++ -template +template void computeHessian( - const std::vector &vec, - std::function(std::vector &)> foo, + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, RowIterator first, RowIterator last, - xad::Tape> *tape = xad::Tape>::getActive()) + Tape> *tape = Tape>::getActive()) ``` This mode uses a [Tape](ref/tape.md) to compute second derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. -#### `fwd_fwd` +### Forward over Forward Mode ```c++ -template +template void computeHessian( - const std::vector &vec, - std::function(std::vector &)> foo, + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, RowIterator first, RowIterator last) ``` This mode does not require a Tape which can help reduce the overhead that comes -with one. +with it, at the expense of requiring more executions of the function given to +determine the full Hessian. ## Example Use -Given $f(x, y, z, w) = sin(x * y) - cos(y * z) - sin(z * w) - cos(w * x)$, or +Given $f(x, y, z, w) = sin(x y) - cos(y z) - sin(z w) - cos(w x)$, or ```c++ auto foo = [](std::vector &x) -> AD @@ -60,7 +62,7 @@ auto foo = [](std::vector &x) -> AD }; ``` -with the derivatives of interest being +with the derivatives of interest calculated at the point ```c++ std::vector x_ad({1.0, 1.5, 1.3, 1.2}); @@ -104,11 +106,11 @@ Note that if no tape is setup, one will be created when computing the Hessian. is define our input values and our function, then call `computeHessian()`: ```c++ - auto foo = [](std::vector &x) -> AD + std::function &)> foo = [](std::vector &x) -> AD { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; - auto hessian = xad::computeHessian(x_ad, foo); + auto hessian = computeHessian(x_ad, foo); ``` Note the signature of `foo()`. Any other signature will throw an error. diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index 9006f9a..d9dea39 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -11,54 +11,60 @@ Jacobians can be computed in `adj` or `fwd` mode. The `computeJacobian()` method takes a set of variables packaged in a `std::vector` and a function with signature -`std::vector foo(std::vector)`. +`std::vector foo(std::vector)`, where `T` is either a forward-mode +or adjoint-mode active type (`FReal` or `AReal`). ## Return Types -If provided with `RowIterators`, `computeHessian()` will write directly to -them and return `void`. If no `RowIterators` are provided, the Hessian will be -written to a `std::vector>` and returned. +If provided with `RowIterators`, `computeJacobian()` will write directly to +them and return `void`. If no `RowIterators` are provided, the Jacobian will be +written to a `std::vector>` and returned, where `T` is the +underlying passive type (usually `double`). ## Specialisations -#### `adj` +### Adjoint Mode ```c++ -template -void computeJacobian(const std::vector &vec, - std::function(std::vector &)> foo, +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, RowIterator first, RowIterator last, - xad::Tape *tape = xad::Tape::getActive()) + Tape *tape = Tape::getActive()) ``` This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will be instantiated within the method or set to the current active Tape using `Tape::getActive()` if none is passed as argument. -#### `fwd_fwd` +### Forward Mode ```c++ -template -void computeJacobian(const std::vector &vec, - std::function(std::vector &)> foo, +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, RowIterator first, RowIterator last) ``` This mode does not require a Tape and can help reduce the overhead that -comes with one. +comes with one. It is recommended for functions that have a higher number +of outputs than inputs. ## Example Use Given $f(x, y, z, w) = [sin(x + y) sin(y + z) cos(z + w) cos(w + x)]$, or ```c++ -auto foo = [](std::vector &x) -> std::vector +std::function(std::vector&)> foo = [](std::vector &x) -> std::vector { - return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; + return {sin(x[0] + x[1]), + sin(x[1] + x[2]), + cos(x[2] + x[3]), + cos(x[3] + x[0])}; }; ``` -with the derivatives of interest being +with the derivatives calculated at the following point ```c++ std::vector x_ad({1.0, 1.5, 1.3, 1.2}); @@ -102,13 +108,13 @@ Note that if no tape is setup, one will be created when computing the Jacobian. define our input values and our function, then call `computeJacobian()`: ```c++ - auto foo = [](std::vector &x) -> std::vector + std::function(std::vector&)> foo = [](std::vector& x) -> std::vector { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; - auto jacobian = xad::computeJacobian(x_ad, foo); + auto jacobian = computeJacobian(x_ad, foo); ``` Note the signature of `foo()`. Any other signature will throw an error. diff --git a/samples/Hessian/main.cpp b/samples/Hessian/main.cpp index c7e1d1e..45aadda 100644 --- a/samples/Hessian/main.cpp +++ b/samples/Hessian/main.cpp @@ -45,10 +45,10 @@ int main() std::vector x_ad({1.0, 1.5, 1.3, 1.2}); // many-input function - auto foo = [](std::vector &x) -> AD + std::function &)> foo = [](std::vector &x) -> AD { return sin(x[0] * x[1]) - cos(x[1] * x[2]) - sin(x[2] * x[3]) - cos(x[3] * x[0]); }; - auto hessian = xad::computeHessian(x_ad, foo); + auto hessian = computeHessian(x_ad, foo); // output results for (auto row : hessian) diff --git a/samples/Jacobian/main.cpp b/samples/Jacobian/main.cpp index df57967..c6e5a7a 100644 --- a/samples/Jacobian/main.cpp +++ b/samples/Jacobian/main.cpp @@ -44,10 +44,10 @@ int main() std::vector x_ad({1.0, 1.5, 1.3, 1.2}); // many-input many-output function - auto foo = [](std::vector &x) -> std::vector + std::function(std::vector&)> foo = [](std::vector& x) -> std::vector { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]), cos(x[3] + x[0])}; }; - auto jacobian = xad::computeJacobian(x_ad, foo); + auto jacobian = computeJacobian(x_ad, foo); // output results for (auto row : jacobian) diff --git a/src/XAD/Hessian.hpp b/src/XAD/Hessian.hpp index 4888622..3018d6d 100644 --- a/src/XAD/Hessian.hpp +++ b/src/XAD/Hessian.hpp @@ -28,6 +28,7 @@ #include #include +#include #include #include @@ -37,9 +38,9 @@ namespace xad // fwd_adj 2d vector template std::vector> computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo, - xad::Tape> *tape = xad::Tape>::getActive()) + const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + Tape> *tape = Tape>::getActive()) { std::vector> matrix(vec.size(), std::vector(vec.size(), 0.0)); computeHessian(vec, foo, begin(matrix), end(matrix), tape); @@ -47,25 +48,23 @@ std::vector> computeHessian( } // fwd_adj iterator -template -void computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo, - RowIterator first, RowIterator last, - xad::Tape> *tape = xad::Tape>::getActive()) +template +void computeHessian(const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + RowIterator first, RowIterator last, + Tape> *tape = Tape>::getActive()) { unsigned int domain(static_cast(vec.size())); if (std::distance(first, last) != domain || std::distance(first->cbegin(), first->cend()) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert( - xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); - std::unique_ptr>> t; + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + std::unique_ptr>> t; if (!tape) { - t = std::unique_ptr>>(new xad::Tape>()); + t = std::unique_ptr>>(new Tape>()); tape = t.get(); } @@ -77,7 +76,7 @@ void computeHessian( { derivative(value(v[i])) = 1.0; tape->newRecording(); - xad::AReal> y = foo(v); + AReal> y = foo(v); tape->registerOutput(y); value(derivative(y)) = 1.0; tape->computeAdjoints(); @@ -94,8 +93,8 @@ void computeHessian( // fwd_fwd 2d vector template std::vector> computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo) + const std::vector>> &vec, + std::function>(std::vector>> &)> foo) { std::vector> matrix(vec.size(), std::vector(vec.size(), 0.0)); computeHessian(vec, foo, begin(matrix), end(matrix)); @@ -103,20 +102,18 @@ std::vector> computeHessian( } // fwd_fwd iterator -template -void computeHessian( - const std::vector>> &vec, - std::function>(std::vector>> &)> foo, - RowIterator first, RowIterator last) +template +void computeHessian(const std::vector>> &vec, + std::function>(std::vector>> &)> foo, + RowIterator first, RowIterator last) { unsigned int domain(static_cast(vec.size())); if (std::distance(first, last) != domain || std::distance(first->cbegin(), first->cend()) != domain) throw OutOfRange("Iterator not allocated enough space"); - static_assert( - xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); auto v(vec); @@ -128,7 +125,7 @@ void computeHessian( for (unsigned int j = 0; j < domain; j++, col++) { derivative(value(v[j])) = 1.0; - xad::FReal> y = foo(v); + FReal> y = foo(v); *col = derivative(derivative(y)); derivative(value(v[j])) = 0.0; } diff --git a/src/XAD/Jacobian.hpp b/src/XAD/Jacobian.hpp index 1e66f3f..1c2778b 100644 --- a/src/XAD/Jacobian.hpp +++ b/src/XAD/Jacobian.hpp @@ -37,9 +37,9 @@ namespace xad // adj 2d vector template std::vector> computeJacobian( - const std::vector> &vec, - std::function>(std::vector> &)> foo, - xad::Tape *tape = xad::Tape::getActive()) + const std::vector> &vec, + std::function>(std::vector> &)> foo, + Tape *tape = Tape::getActive()) { auto v(vec); std::vector> matrix(foo(v).size(), std::vector(v.size(), 0.0)); @@ -48,21 +48,19 @@ std::vector> computeJacobian( } // adj iterator -template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, - RowIterator first, RowIterator last, - xad::Tape *tape = xad::Tape::getActive()) +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last, Tape *tape = Tape::getActive()) { if (std::distance(first->cbegin(), first->cend()) != vec.size()) throw OutOfRange("Iterator not allocated enough space (domain)"); - static_assert( - xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); - std::unique_ptr> t; + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); + std::unique_ptr> t; if (!tape) { - t = std::unique_ptr>(new xad::Tape()); + t = std::unique_ptr>(new Tape()); tape = t.get(); } @@ -91,8 +89,8 @@ void computeJacobian(const std::vector> &vec, // fwd 2d vector template std::vector> computeJacobian( - const std::vector> &vec, - std::function>(std::vector> &)> foo) + const std::vector> &vec, + std::function>(std::vector> &)> foo) { auto v(vec); std::vector> matrix(foo(v).size(), std::vector(v.size(), 0.0)); @@ -101,16 +99,15 @@ std::vector> computeJacobian( } // fwd iterator -template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, +template +void computeJacobian(const std::vector> &vec, + std::function>(std::vector> &)> foo, RowIterator first, RowIterator last) { if (std::distance(first->cbegin(), first->cend()) != vec.size()) throw OutOfRange("Iterator not allocated enough space (domain)"); - static_assert( - xad::detail::has_begin::value_type>::value, - "RowIterator must dereference to a type that implements a begin() method"); + static_assert(detail::has_begin::value_type>::value, + "RowIterator must dereference to a type that implements a begin() method"); auto v(vec); unsigned int domain = static_cast(vec.size()), diff --git a/test/Hessian_test.cpp b/test/Hessian_test.cpp index 0749db7..61a444e 100644 --- a/test/Hessian_test.cpp +++ b/test/Hessian_test.cpp @@ -49,13 +49,12 @@ TEST(HessianTest, QuadraticForwardAdjoint) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticForwardAdjointAutoTape) @@ -68,13 +67,12 @@ TEST(HessianTest, QuadraticForwardAdjointAutoTape) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo); + std::vector> computed_hessian = xad::computeHessian(x, foo); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticForwardAdjointFetchTape) @@ -90,13 +88,12 @@ TEST(HessianTest, QuadraticForwardAdjointFetchTape) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::vector> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo); + std::vector> computed_hessian = xad::computeHessian(x, foo); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticForwardAdjointWithIterator) @@ -112,25 +109,17 @@ TEST(HessianTest, QuadraticForwardAdjointWithIterator) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::list> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; + std::list> expected_hessian = {{2.0, 0.0}, {0.0, 2.0}}; - std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); xad::computeHessian(x, foo, begin(computed_hessian), end(computed_hessian), &tape); - auto row1 = computed_hessian.begin(), row2 = expected_hessian.begin(); - while (row1 != computed_hessian.end() && row2 != expected_hessian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_hessian.begin(); + auto computed_it = computed_hessian.begin(); + for (; expected_it != expected_hessian.end() && computed_it != computed_hessian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(HessianTest, SingleInputForwardAdjoint) @@ -146,13 +135,12 @@ TEST(HessianTest, SingleInputForwardAdjoint) // f(x) = x[0]^3 + x[0] auto foo = [](std::vector &x) -> AD { return x[0] * x[0] * x[0] + x[0]; }; - std::vector> expected_hessian = {{18.0}}; + std::vector> expected_hessian = {{18.0}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticForwardForward) @@ -165,14 +153,13 @@ TEST(HessianTest, QuadraticForwardForward) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::vector> expected_hessian = {{2.0, 0.0}, // - {0.0, 2.0}}; + std::vector> expected_hessian = {{2.0, 0.0}, // + {0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo); + std::vector> computed_hessian = xad::computeHessian(x, foo); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticForwardForwardWithIterator) @@ -185,26 +172,18 @@ TEST(HessianTest, QuadraticForwardForwardWithIterator) // f(x) = x[0]^2 + x[1]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1]; }; - std::list> expected_hessian = {{2.0, 0.0}, // - {0.0, 2.0}}; + std::list> expected_hessian = {{2.0, 0.0}, // + {0.0, 2.0}}; - std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); + std::list> computed_hessian(x.size(), std::list(x.size(), 0.0)); xad::computeHessian(x, foo, begin(computed_hessian), end(computed_hessian)); - auto row1 = computed_hessian.begin(), row2 = expected_hessian.begin(); - while (row1 != computed_hessian.end() && row2 != expected_hessian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_hessian.begin(); + auto computed_it = computed_hessian.begin(); + for (; expected_it != expected_hessian.end() && computed_it != computed_hessian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(HessianTest, SingleInputForwardForward) @@ -217,13 +196,12 @@ TEST(HessianTest, SingleInputForwardForward) // f(x) = x[0]^3 + x[0] auto foo = [](std::vector &x) -> AD { return x[0] * x[0] * x[0] + x[0]; }; - std::vector> expected_hessian = {{18.0}}; + std::vector> expected_hessian = {{18.0}}; - auto computed_hessian = xad::computeHessian(x, foo); + std::vector> computed_hessian = xad::computeHessian(x, foo); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticThreeVariablesForwardAdjoint) @@ -239,14 +217,13 @@ TEST(HessianTest, QuadraticThreeVariablesForwardAdjoint) // f(x) = x[0]^2 + x[1]^2 + x[2]^2 auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1] + x[2] * x[2]; }; - std::vector> expected_hessian = { + std::vector> expected_hessian = { {2.0, 0.0, 0.0}, {0.0, 2.0, 0.0}, {0.0, 0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, ComplexFunctionForwardAdjoint) @@ -262,16 +239,16 @@ TEST(HessianTest, ComplexFunctionForwardAdjoint) // f(x) = x[0] * sin(x[1]) + x[2] * exp(x[3]) auto foo = [](std::vector &x) -> AD { return x[0] * sin(x[1]) + x[2] * exp(x[3]); }; - std::vector> expected_hessian = {{0.0, cos(x[1]), 0.0, 0.0}, - {cos(x[1]), -x[0] * sin(x[1]), 0.0, 0.0}, - {0.0, 0.0, 0.0, exp(x[3])}, - {0.0, 0.0, exp(x[3]), x[2] * exp(x[3])}}; + std::vector> expected_hessian = { + {0.0, cos(value(value(x[1]))), 0.0, 0.0}, + {cos(value(value(x[1]))), -value(value(x[0])) * sin(value(value(x[1]))), 0.0, 0.0}, + {0.0, 0.0, 0.0, exp(value(value(x[3])))}, + {0.0, 0.0, exp(value(value(x[3]))), value(value(x[2])) * exp(value(value(x[3])))}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, FourthOrderPolynomialForwardAdjoint) @@ -287,15 +264,15 @@ TEST(HessianTest, FourthOrderPolynomialForwardAdjoint) // f(x) = x[0]^4 + x[1]^4 + x[2]^4 auto foo = [](std::vector &x) -> AD { return pow(x[0], 4) + pow(x[1], 4) + pow(x[2], 4); }; - std::vector> expected_hessian = {{12.0 * x[0] * x[0], 0.0, 0.0}, - {0.0, 12.0 * x[1] * x[1], 0.0}, - {0.0, 0.0, 12.0 * x[2] * x[2]}}; + std::vector> expected_hessian = { + {12.0 * value(value(x[0])) * value(value(x[0])), 0.0, 0.0}, + {0.0, 12.0 * value(value(x[1])) * value(value(x[1])), 0.0}, + {0.0, 0.0, 12.0 * value(value(x[2])) * value(value(x[2]))}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, HigherOrderInteractionForwardAdjoint) @@ -311,15 +288,15 @@ TEST(HessianTest, HigherOrderInteractionForwardAdjoint) // f(x) = x[0] * x[1] * x[2] auto foo = [](std::vector &x) -> AD { return x[0] * x[1] * x[2]; }; - std::vector> expected_hessian = {{0.0, x[2], x[1]}, // - {x[2], 0.0, x[0]}, - {x[1], x[0], 0.0}}; + std::vector> expected_hessian = { + {0.0, value(value(x[2])), value(value(x[1]))}, // + {value(value(x[2])), 0.0, value(value(x[0]))}, + {value(value(x[1])), value(value(x[0])), 0.0}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, QuadraticFourVariablesForwardAdjoint) @@ -336,16 +313,15 @@ TEST(HessianTest, QuadraticFourVariablesForwardAdjoint) auto foo = [](std::vector &x) -> AD { return x[0] * x[0] + x[1] * x[1] + x[2] * x[2] + x[3] * x[3]; }; - std::vector> expected_hessian = {{2.0, 0.0, 0.0, 0.0}, // - {0.0, 2.0, 0.0, 0.0}, - {0.0, 0.0, 2.0, 0.0}, - {0.0, 0.0, 0.0, 2.0}}; + std::vector> expected_hessian = {{2.0, 0.0, 0.0, 0.0}, // + {0.0, 2.0, 0.0, 0.0}, + {0.0, 0.0, 2.0, 0.0}, + {0.0, 0.0, 0.0, 2.0}}; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, LargeHessianForwardAdjoint) @@ -369,14 +345,13 @@ TEST(HessianTest, LargeHessianForwardAdjoint) return result; }; - std::vector> expected_hessian(16, std::vector(16, 1.0)); + std::vector> expected_hessian(16, std::vector(16, 1.0)); for (size_t i = 0; i < 16; ++i) expected_hessian[i][i] = 2.0; - auto computed_hessian = xad::computeHessian(x, foo, &tape); + std::vector> computed_hessian = xad::computeHessian(x, foo, &tape); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, LargeHessianForwardForward) @@ -397,14 +372,13 @@ TEST(HessianTest, LargeHessianForwardForward) return result; }; - std::vector> expected_hessian(16, std::vector(16, 1.0)); + std::vector> expected_hessian(16, std::vector(16, 1.0)); for (size_t i = 0; i < 16; ++i) expected_hessian[i][i] = 2.0; - auto computed_hessian = xad::computeHessian(x, foo); + std::vector> computed_hessian = xad::computeHessian(x, foo); for (unsigned int i = 0; i < expected_hessian.size(); i++) - for (unsigned int j = 0; j < expected_hessian[i].size(); j++) - ASSERT_EQ(expected_hessian[i][j], computed_hessian[i][j]); + EXPECT_THAT(computed_hessian[i], Pointwise(DoubleEq(), expected_hessian[i])); } TEST(HessianTest, OutOfBoundsDomainSizeMismatch) @@ -416,14 +390,14 @@ TEST(HessianTest, OutOfBoundsDomainSizeMismatch) auto foo = [](std::vector &x) -> AD { return x[0]; }; - std::vector> jacobian(2, std::vector(3)); + std::vector> jacobian(2, std::vector(3)); auto launch = [](std::vector x, std::function>( std::vector>> &)> foo, - std::vector>>>::iterator first, - std::vector>>>::iterator last) + std::vector>::iterator first, + std::vector>::iterator last) { using RowIterator = decltype(first); xad::computeHessian(x, foo, first, last); diff --git a/test/Jacobian_test.cpp b/test/Jacobian_test.cpp index 12bcca2..bdc958c 100644 --- a/test/Jacobian_test.cpp +++ b/test/Jacobian_test.cpp @@ -33,6 +33,8 @@ #include #include +using namespace ::testing; + TEST(JacobianTest, SimpleAdjoint) { typedef xad::adj mode; @@ -47,14 +49,14 @@ TEST(JacobianTest, SimpleAdjoint) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::vector> expected_jacobian = {{1.0, cos(x[1])}, // - {cos(x[0]), 1.0}}; + std::vector> expected_jacobian = {{1.0, cos(value(x[1]))}, // + {cos(value(x[0])), 1.0}}; - auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + std::vector> computed_jacobian = + xad::computeJacobian(x, foo, &tape); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, SimpleAdjointIteratorAutoTape) @@ -68,26 +70,18 @@ TEST(JacobianTest, SimpleAdjointIteratorAutoTape) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::list> expected_jacobian = {{1.0, cos(x[1])}, // - {cos(x[0]), 1.0}}; + std::list> expected_jacobian = {{1.0, cos(value(x[1]))}, // + {cos(value(x[0])), 1.0}}; - std::list> computed_jacobian(2, std::list(2, 0.0)); + std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_jacobian.begin(); + auto computed_it = computed_jacobian.begin(); + for (; expected_it != expected_jacobian.end() && computed_it != computed_jacobian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(JacobianTest, SimpleAdjointIteratorFetchTape) @@ -104,26 +98,18 @@ TEST(JacobianTest, SimpleAdjointIteratorFetchTape) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::list> expected_jacobian = {{1.0, cos(x[1])}, // - {cos(x[0]), 1.0}}; + std::list> expected_jacobian = {{1.0, cos(value(x[1]))}, // + {cos(value(x[0])), 1.0}}; - std::list> computed_jacobian(2, std::list(2, 0.0)); + std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_jacobian.begin(); + auto computed_it = computed_jacobian.begin(); + for (; expected_it != expected_jacobian.end() && computed_it != computed_jacobian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(JacobianTest, SimpleForward) @@ -137,14 +123,13 @@ TEST(JacobianTest, SimpleForward) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::vector> expected_jacobian = {{1.0, cos(x[1])}, // - {cos(x[0]), 1.0}}; + std::vector> expected_jacobian = {{1.0, cos(value(x[1]))}, // + {cos(value(x[0])), 1.0}}; - auto computed_jacobian = xad::computeJacobian(x, foo); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, SimpleForwardIterator) @@ -158,26 +143,18 @@ TEST(JacobianTest, SimpleForwardIterator) auto foo = [](std::vector &x) -> std::vector { return {x[0] + sin(x[1]), x[1] + sin(x[0])}; }; - std::list> expected_jacobian = {{1.0, cos(x[1])}, // - {cos(x[0]), 1.0}}; + std::list> expected_jacobian = {{1.0, cos(value(x[1]))}, // + {cos(value(x[0])), 1.0}}; - std::list> computed_jacobian(2, std::list(2, 0.0)); + std::list> computed_jacobian(2, std::list(2, 0.0)); xad::computeJacobian( x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_jacobian.begin(); + auto computed_it = computed_jacobian.begin(); + for (; expected_it != expected_jacobian.end() && computed_it != computed_jacobian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(JacobianTest, ComplexFunctionAdjoint) @@ -194,14 +171,15 @@ TEST(JacobianTest, ComplexFunctionAdjoint) auto foo = [](std::vector &x) -> std::vector { return {x[0] * x[1], x[2] * exp(x[3])}; }; - std::vector> expected_jacobian = {{x[1], x[0], 0.0, 0.0}, - {0.0, 0.0, exp(x[3]), x[2] * exp(x[3])}}; + std::vector> expected_jacobian = { + {value(x[1]), value(x[0]), 0.0, 0.0}, + {0.0, 0.0, exp(value(x[3])), value(x[2]) * exp(value(x[3]))}}; - auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + std::vector> computed_jacobian = + xad::computeJacobian(x, foo, &tape); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, DomainLargerThanCodomainForward) @@ -214,14 +192,13 @@ TEST(JacobianTest, DomainLargerThanCodomainForward) // f(x) = [ x[0] + x[1], x[2] * x[3] ] auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[2] * x[3]}; }; - std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // - {0.0, 0.0, x[3], x[2]}}; + std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // + {0.0, 0.0, value(x[3]), value(x[2])}}; - auto computed_jacobian = xad::computeJacobian(x, foo); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, DomainLargerThanCodomainAdjoint) @@ -234,14 +211,13 @@ TEST(JacobianTest, DomainLargerThanCodomainAdjoint) // f(x) = [ x[0] + x[1], x[2] * x[3] ] auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[2] * x[3]}; }; - std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // - {0.0, 0.0, x[3], x[2]}}; + std::vector> expected_jacobian = {{1.0, 1.0, 0.0, 0.0}, // + {0.0, 0.0, value(x[3]), value(x[2])}}; - auto computed_jacobian = xad::computeJacobian(x, foo); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, DomainSmallerThanCodomainAdjoint) @@ -258,15 +234,15 @@ TEST(JacobianTest, DomainSmallerThanCodomainAdjoint) auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[0] - x[1], x[0] * x[1]}; }; - std::vector> expected_jacobian = {{1.0, 1.0}, // - {1.0, -1.0}, - {x[1], x[0]}}; + std::vector> expected_jacobian = {{1.0, 1.0}, // + {1.0, -1.0}, + {value(x[1]), value(x[0])}}; - auto computed_jacobian = xad::computeJacobian(x, foo, &tape); + std::vector> computed_jacobian = + xad::computeJacobian(x, foo, &tape); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, ComplexDomainNotEqualCodomainForwardIterator) @@ -280,27 +256,19 @@ TEST(JacobianTest, ComplexDomainNotEqualCodomainForwardIterator) auto foo = [](std::vector &x) -> std::vector { return {x[0] + x[1], x[1] * x[2], exp(x[0])}; }; - std::list> expected_jacobian = {{1.0, 1.0, 0.0}, // - {0.0, x[2], x[1]}, - {exp(x[0]), 0.0, 0.0}}; + std::list> expected_jacobian = {{1.0, 1.0, 0.0}, // + {0.0, value(x[2]), value(x[1])}, + {exp(value(x[0])), 0.0, 0.0}}; - std::list> computed_jacobian(3, std::list(3, 0.0)); + std::list> computed_jacobian(3, std::list(3, 0.0)); xad::computeJacobian( x, foo, begin(computed_jacobian), end(computed_jacobian)); - auto row1 = computed_jacobian.begin(), row2 = expected_jacobian.begin(); - while (row1 != computed_jacobian.end() && row2 != expected_jacobian.end()) - { - auto col1 = row1->begin(), col2 = row2->begin(); - while (col1 != row1->end() && col2 != row2->end()) - { - ASSERT_EQ(*col1, *col2); - col1++; - col2++; - } - row1++; - row2++; - } + auto expected_it = expected_jacobian.begin(); + auto computed_it = computed_jacobian.begin(); + for (; expected_it != expected_jacobian.end() && computed_it != computed_jacobian.end(); + ++expected_it, ++computed_it) + EXPECT_THAT(*computed_it, Pointwise(DoubleEq(), *expected_it)); } TEST(JacobianTest, TrigonometricFunctionForward) @@ -313,14 +281,13 @@ TEST(JacobianTest, TrigonometricFunctionForward) // f(x) = [ sin(x[0]), cos(x[1]) ] auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; - std::vector> expected_jacobian = {{cos(x[0]), 0.0}, // - {0.0, -sin(x[1])}}; + std::vector> expected_jacobian = {{cos(value(x[0])), 0.0}, // + {0.0, -sin(value(x[1]))}}; - auto computed_jacobian = xad::computeJacobian(x, foo); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, TrigonometricFunctionAdjointAutoTape) @@ -333,16 +300,13 @@ TEST(JacobianTest, TrigonometricFunctionAdjointAutoTape) // f(x) = [ sin(x[0]), cos(x[1]) ] auto foo = [](std::vector &x) -> std::vector { return {sin(x[0]), cos(x[1])}; }; - std::vector> expected_jacobian = {{cos(x[0]), 0.0}, // - {0.0, -sin(x[1])}}; + std::vector> expected_jacobian = {{cos(value(x[0])), 0.0}, // + {0.0, -sin(value(x[1]))}}; - auto computed_jacobian = xad::computeJacobian(x, foo); + std::vector> computed_jacobian = xad::computeJacobian(x, foo); for (unsigned int i = 0; i < expected_jacobian.size(); i++) - for (unsigned int j = 0; j < expected_jacobian[i].size(); j++) - { - ASSERT_EQ(expected_jacobian[i][j], computed_jacobian[i][j]); - } + EXPECT_THAT(computed_jacobian[i], Pointwise(DoubleEq(), expected_jacobian[i])); } TEST(JacobianTest, OutOfBoundsDomainSizeMismatch) @@ -354,13 +318,13 @@ TEST(JacobianTest, OutOfBoundsDomainSizeMismatch) auto foo = [](std::vector &x) -> std::vector { return {x[0], x[1]}; }; - std::vector> jacobian(2, std::vector(3)); + std::vector> jacobian(2, std::vector(3)); auto launch = [](std::vector x, std::function>(std::vector> &)> foo, - std::vector>>::iterator first, - std::vector>>::iterator last) + std::vector>::iterator first, + std::vector>::iterator last) { using RowIterator = decltype(first); xad::computeJacobian(x, foo, first, last); @@ -378,13 +342,13 @@ TEST(JacobianTest, OutOfBoundsCodomainSizeMismatch) auto foo = [](std::vector &x) -> std::vector { return {x[0], x[0]}; }; - std::vector> jacobian(1, std::vector(1)); + std::vector> jacobian(1, std::vector(1)); auto launch = [](std::vector x, std::function>(std::vector> &)> foo, - std::vector>>::iterator first, - std::vector>>::iterator last) + std::vector>::iterator first, + std::vector>::iterator last) { using RowIterator = decltype(first); xad::computeJacobian(x, foo, first, last); From a6263b670a7448dbbe2a2284882f9232e94afa2d Mon Sep 17 00:00:00 2001 From: imraneamri Date: Wed, 12 Jun 2024 10:49:36 +0100 Subject: [PATCH 30/30] Final changes for Codacy (the rest require config) --- docs/ref/jacobian.md | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/docs/ref/jacobian.md b/docs/ref/jacobian.md index d9dea39..13f47ff 100644 --- a/docs/ref/jacobian.md +++ b/docs/ref/jacobian.md @@ -27,10 +27,11 @@ underlying passive type (usually `double`). ```c++ template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, - RowIterator first, RowIterator last, - Tape *tape = Tape::getActive()) +void computeJacobian( + const std::vector> &vec, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last, + Tape *tape = Tape::getActive()) ``` This mode uses a [Tape](ref/tape.md) to compute derivatives. This Tape will @@ -41,9 +42,10 @@ be instantiated within the method or set to the current active Tape using ```c++ template -void computeJacobian(const std::vector> &vec, - std::function>(std::vector> &)> foo, - RowIterator first, RowIterator last) +void computeJacobian( + const std::vector> &vec, + std::function>(std::vector> &)> foo, + RowIterator first, RowIterator last) ``` This mode does not require a Tape and can help reduce the overhead that @@ -55,8 +57,8 @@ of outputs than inputs. Given $f(x, y, z, w) = [sin(x + y) sin(y + z) cos(z + w) cos(w + x)]$, or ```c++ -std::function(std::vector&)> foo = [](std::vector &x) -> std::vector -{ +std::function(std::vector&)> foo = +[](std::vector &x) -> std::vector { return {sin(x[0] + x[1]), sin(x[1] + x[2]), cos(x[2] + x[3]),