From 650e5cc6b827d3d51fc5a1c9c935f7b650e53c08 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Valentin=20Leon=20St=C3=B6cker?= Date: Tue, 20 Jan 2026 17:30:34 +0100 Subject: [PATCH] init --- .gitignore | 9 + CMakeLists.txt | 20 ++ flake.lock | 61 ++++++ flake.nix | 71 +++++++ include/autoopt/derivative.hpp | 61 ++++++ include/autoopt/dual.hpp | 241 +++++++++++++++++++++++ include/autoopt/ellipse.hpp | 50 +++++ include/autoopt/optimization_problem.hpp | 111 +++++++++++ include/autoopt/quadric.hpp | 60 ++++++ include/autoopt/util.hpp | 38 ++++ tests/CMakeLists.txt | 10 + tests/dual.cpp | 137 +++++++++++++ tests/ellipse.cpp | 78 ++++++++ tests/quadric.cpp | 21 ++ 14 files changed, 968 insertions(+) create mode 100644 .gitignore create mode 100644 CMakeLists.txt create mode 100644 flake.lock create mode 100644 flake.nix create mode 100644 include/autoopt/derivative.hpp create mode 100644 include/autoopt/dual.hpp create mode 100644 include/autoopt/ellipse.hpp create mode 100644 include/autoopt/optimization_problem.hpp create mode 100644 include/autoopt/quadric.hpp create mode 100644 include/autoopt/util.hpp create mode 100644 tests/CMakeLists.txt create mode 100644 tests/dual.cpp create mode 100644 tests/ellipse.cpp create mode 100644 tests/quadric.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..eb2c859 --- /dev/null +++ b/.gitignore @@ -0,0 +1,9 @@ +# nix +result + +# cmake +build + +# direnv files +.direnv +.envrc \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..116105a --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,20 @@ +cmake_minimum_required(VERSION 4.0) + +project(autoopt) + +set(CMAKE_CXX_STANDARD 23) + +add_compile_options(-Wall -Werror -Wpedantic) + +add_library(autoopt INTERFACE) +target_include_directories(autoopt INTERFACE include) + +install(DIRECTORY include/autoopt DESTINATION include) + +# add tests if gtest is found +find_library(GTestPackage gtest QUIET) +if(GTestPackage) + enable_testing() + add_subdirectory(tests) +endif() + diff --git a/flake.lock b/flake.lock new file mode 100644 index 0000000..5a799b0 --- /dev/null +++ b/flake.lock @@ -0,0 +1,61 @@ +{ + "nodes": { + "flake-utils": { + "inputs": { + "systems": "systems" + }, + "locked": { + "lastModified": 1731533236, + "narHash": "sha256-l0KFg5HjrsfsO/JpG+r7fRrqm12kzFHyUHqHCVpMMbI=", + "owner": "numtide", + "repo": "flake-utils", + "rev": "11707dc2f618dd54ca8739b309ec4fc024de578b", + "type": "github" + }, + "original": { + "owner": "numtide", + "repo": "flake-utils", + "type": "github" + } + }, + "nixpkgs": { + "locked": { + "lastModified": 1768564909, + "narHash": "sha256-Kell/SpJYVkHWMvnhqJz/8DqQg2b6PguxVWOuadbHCc=", + "owner": "nixos", + "repo": "nixpkgs", + "rev": "e4bae1bd10c9c57b2cf517953ab70060a828ee6f", + "type": "github" + }, + "original": { + "owner": "nixos", + "ref": "nixos-unstable", + "repo": "nixpkgs", + "type": "github" + } + }, + "root": { + "inputs": { + "flake-utils": "flake-utils", + "nixpkgs": "nixpkgs" + } + }, + "systems": { + "locked": { + "lastModified": 1681028828, + "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", + "owner": "nix-systems", + "repo": "default", + "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", + "type": "github" + }, + "original": { + "owner": "nix-systems", + "repo": "default", + "type": "github" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix new file mode 100644 index 0000000..110641d --- /dev/null +++ b/flake.nix @@ -0,0 +1,71 @@ +{ + description = "auto-opt-cpp"; + + inputs = { + nixpkgs.url = "github:nixos/nixpkgs?ref=nixos-unstable"; + flake-utils.url = "github:numtide/flake-utils"; + }; + + outputs = + { + self, + nixpkgs, + flake-utils, + }: + flake-utils.lib.eachDefaultSystem ( + system: + let + pkgs = nixpkgs.legacyPackages.${system}; + pkg = + { + stdenv, + cmake, + ninja, + lib, + gtest, + withTests ? false, + }: + let + fs = lib.fileset; + in + stdenv.mkDerivation { + pname = "autoopt-cpp"; + version = "0.1.0"; + + src = fs.toSource { + root = ./.; + fileset = fs.unions [ + ./include + ./CMakeLists.txt + ./tests + ]; + }; + + buildInputs = lib.optionals withTests [ gtest ]; + + nativeBuildInputs = [ + cmake + ninja + ]; + }; + in + rec { + packages = rec { + autoopt-cpp = pkgs.callPackage pkg { }; + autoopt-cpp-with-tests = pkgs.callPackage pkg { withTests = true; }; + default = autoopt-cpp; + }; + devShells = { + default = pkgs.mkShell { + inputsFrom = [ + packages.autoopt-cpp-with-tests + ]; + + shellHook = '' + echo "Development shell for autoopt-cpp" + ''; + }; + }; + } + ); +} diff --git a/include/autoopt/derivative.hpp b/include/autoopt/derivative.hpp new file mode 100644 index 0000000..73ea64d --- /dev/null +++ b/include/autoopt/derivative.hpp @@ -0,0 +1,61 @@ +#pragma once + +#include + +#include "autoopt/dual.hpp" + +namespace autoopt { + +template +T derivative(Func&& f, const T& x) { + dual a(x, T(1)); + dual b = f(a); + return b._dx; +} + +template +std::array gradient(Func&& f, const std::array& x) { + std::array grad{}; + std::array, N> dual_x{}; + for (std::size_t i = 0; i < N; ++i) { + dual_x[i] = dual(x[i], T(0)); + } + for (std::size_t i = 0; i < N; ++i) { + dual_x[i]._dx = T(1); + dual dual_y = f(dual_x); + grad[i] = dual_y._dx; + dual_x[i]._dx = T(0); + } + return grad; +} + +template +using matrix_t = std::array, N>; + +template +matrix_t jacobian(Func&& f, const std::array& x) { + matrix_t jacob{}; + std::array, N> dual_x{}; + for (std::size_t i = 0; i < N; ++i) { + dual_x[i] = dual(x[i], T(0)); + } + for (std::size_t i = 0; i < N; ++i) { + dual_x[i]._dx = T(1); + std::array, M> dual_y = f(dual_x); + for (std::size_t j = 0; j < M; ++j) { + jacob[j][i] = dual_y[j]._dx; + } + dual_x[i]._dx = T(0); + } + return jacob; +} + +template +matrix_t hessian(Func&& f, const std::array& x) { + auto helper_func = [&f](const std::array& y) { + return gradient(f, y); + }; + return jacobian(helper_func, x); +} + +} // namespace autoopt \ No newline at end of file diff --git a/include/autoopt/dual.hpp b/include/autoopt/dual.hpp new file mode 100644 index 0000000..11839b8 --- /dev/null +++ b/include/autoopt/dual.hpp @@ -0,0 +1,241 @@ +#pragma once + +#include + +namespace autoopt { + +template +struct dual { + T _x; + T _dx; + + constexpr dual(T x = T{0}, T dx = T{0}) : _x(x), _dx(dx) {} + + // allow arbitrary deeply-nested dual construction + template + constexpr dual(U x) : _x(T(x)), _dx(T(0)) {} +}; + +template +constexpr dual operator+(const dual& a, const dual& b) { + return dual(a._x + b._x, a._dx + b._dx); +} + +template +constexpr dual operator-(const dual& a) { + return dual(-a._x, -a._dx); +} + +template +constexpr dual operator-(const dual& a, const dual& b) { + return dual(a._x - b._x, a._dx - b._dx); +} + +template +constexpr dual operator*(const dual& a, const dual& b) { + return dual(a._x * b._x, a._x * b._dx + a._dx * b._x); +} + +template +constexpr dual operator/(const dual& a, const dual& b) { + return dual(a._x / b._x, (a._dx * b._x - a._x * b._dx) / (b._x * b._x)); +} + +template +constexpr bool operator==(const dual& a, const dual& b) { + return a._x == b._x; +} + +template +constexpr bool operator!=(const dual& a, const dual& b) { + return a._x != b._x; +} + +template +constexpr bool operator<(const dual& a, const dual& b) { + return a._x < b._x; +} + +template +constexpr bool operator<=(const dual& a, const dual& b) { + return a._x <= b._x; +} + +template +constexpr bool operator>(const dual& a, const dual& b) { + return a._x > b._x; +} + +template +constexpr bool operator>=(const dual& a, const dual& b) { + return a._x >= b._x; +} + +} // namespace autoopt + +namespace std { + +using autoopt::dual; + +// forward declarations of standard functions + +template +constexpr dual abs(const dual& a); + +template +constexpr dual exp(const dual& a); + +template +constexpr dual log(const dual& a); + +template +constexpr dual pow(const dual& a, const dual& b); + +template +constexpr dual sqrt(const dual& a); + +// forward declarations of trigonometric functions + +template +constexpr dual sin(const dual& a); + +template +constexpr dual cos(const dual& a); + +template +constexpr dual tan(const dual& a); + +template +constexpr dual asin(const dual& a); + +template +constexpr dual acos(const dual& a); + +template +constexpr dual atan(const dual& a); + +template +constexpr dual atan2(const dual& y, const dual& x); + +// forward declarations of hyperbolic functions + +template +constexpr dual sinh(const dual& a); + +template +constexpr dual cosh(const dual& a); + +template +constexpr dual tanh(const dual& a); + +template +constexpr dual asinh(const dual& a); + +template +constexpr dual acosh(const dual& a); + +template +constexpr dual atanh(const dual& a); + +// standard functions + +template +constexpr dual abs(const dual& a) { + return dual(std::abs(a._x), (a._x >= T(0) ? T(1) : T(-1)) * a._dx); +} + +template +constexpr dual exp(const dual& a) { + T exp_x = std::exp(a._x); + return dual(exp_x, exp_x * a._dx); +} + +template +constexpr dual log(const dual& a) { + return dual(std::log(a._x), (T(1) / a._x) * a._dx); +} + +template +constexpr dual pow(const dual& a, const dual& b) { + return std::exp(b * std::log(a)); +} + +template +constexpr dual sqrt(const dual& a) { + T sqrt_x = std::sqrt(a._x); + return dual(sqrt_x, (T(1) / (T(2) * sqrt_x)) * a._dx); +} + +// trigonometric functions + +template +constexpr dual sin(const dual& a) { + return dual(std::sin(a._x), std::cos(a._x) * a._dx); +} + +template +constexpr dual cos(const dual& a) { + return dual(std::cos(a._x), -std::sin(a._x) * a._dx); +} + +template +constexpr dual tan(const dual& a) { + T cos_x = std::cos(a._x); + return dual(std::tan(a._x), (T(1) / (cos_x * cos_x)) * a._dx); +} + +template +constexpr dual asin(const dual& a) { + return dual(std::asin(a._x), (T(1) / std::sqrt(T(1) - a._x * a._x)) * a._dx); +} + +template +constexpr dual acos(const dual& a) { + return dual(std::acos(a._x), (-T(1) / std::sqrt(T(1) - a._x * a._x)) * a._dx); +} + +template +constexpr dual atan(const dual& a) { + return dual(std::atan(a._x), (T(1) / (T(1) + a._x * a._x)) * a._dx); +} + +template +constexpr dual atan2(const dual& y, const dual& x) { + return dual(std::atan2(y._x, x._x), + (x._x * y._dx - y._x * x._dx) / (x._x * x._x + y._x * y._x)); +} + +// hyperbolic functions + +template +constexpr dual sinh(const dual& a) { + return dual(std::sinh(a._x), std::cosh(a._x) * a._dx); +} + +template +constexpr dual cosh(const dual& a) { + return dual(std::cosh(a._x), std::sinh(a._x) * a._dx); +} + +template +constexpr dual tanh(const dual& a) { + T cosh_x = std::cosh(a._x); + return dual(std::tanh(a._x), (T(1) / (cosh_x * cosh_x)) * a._dx); +} + +template +constexpr dual asinh(const dual& a) { + return dual(std::asinh(a._x), (T(1) / std::sqrt(a._x * a._x + T(1))) * a._dx); +} + +template +constexpr dual acosh(const dual& a) { + return dual(std::acosh(a._x), (T(1) / std::sqrt(a._x * a._x - T(1))) * a._dx); +} + +template +constexpr dual atanh(const dual& a) { + return dual(std::atanh(a._x), (T(1) / (T(1) - a._x * a._x)) * a._dx); +} + +} // namespace std diff --git a/include/autoopt/ellipse.hpp b/include/autoopt/ellipse.hpp new file mode 100644 index 0000000..40f49aa --- /dev/null +++ b/include/autoopt/ellipse.hpp @@ -0,0 +1,50 @@ +#pragma once + +#include "autoopt/quadric.hpp" + +namespace autoopt { + +template +struct ellipse { + T left_arm; + T right_arm; + T entrance_angle; + + constexpr quadric to_quadric() const { + T a = (left_arm + right_arm) / T{2}; + T a2 = a * a; + + T left_x = -left_arm * std::cos(entrance_angle); + T left_y = left_arm * std::sin(entrance_angle); + + T right_x = right_arm * std::cos(entrance_angle); + T right_y = right_arm * std::sin(entrance_angle); + + T c_x = (left_x + right_x) / T{2}; + T c_y = (left_y + right_y) / T{2}; + + T c2 = (left_x - c_x) * (left_x - c_x) + (left_y - c_y) * (left_y - c_y); + T b2 = a2 - c2; + + + // source: https://en.wikipedia.org/wiki/Ellipse#General_ellipse + + T theta = std::atan2(right_y - left_y, right_x - left_x); + T cos_t = std::cos(theta); + T sin_t = std::sin(theta); + T cos_t2 = cos_t * cos_t; + T sin_t2 = sin_t * sin_t; + T sico_t = sin_t * cos_t; + + T A = a2 * sin_t2 + b2 * cos_t2; + T B = T{2} * (b2 - a2) * sico_t; + T C = a2 * cos_t2 + b2 * sin_t2; + T D = -T{2} * A * c_x - B * c_y; + T E = -B * c_x - T{2} * C * c_y; + T F = A * c_x * c_x + B * c_x * c_y + C * c_y * c_y - a2 * b2; + + return quadric(A, B, C, D, E, F); + } +}; + +} // namespace autoopt \ No newline at end of file diff --git a/include/autoopt/optimization_problem.hpp b/include/autoopt/optimization_problem.hpp new file mode 100644 index 0000000..8ac78c7 --- /dev/null +++ b/include/autoopt/optimization_problem.hpp @@ -0,0 +1,111 @@ +#pragma once + +#include +#include + +#include "autoopt/derivative.hpp" + +namespace autoopt { + +template +struct optimization_problem { + virtual std::array initial_guess() = 0; + + virtual T objective(const std::array& params) = 0; + virtual std::array gradient(const std::array& params) = 0; + virtual matrix_t hessian(const std::array& params) = 0; +}; + +template +struct auto_diff_optimization_problem : public optimization_problem { + Func _objective_func; + std::array _initial_guess; + + auto_diff_optimization_problem( + Func objective_func, std::array initial_guess = std::array{}) + : _objective_func(objective_func), _initial_guess(initial_guess) {} + + std::array initial_guess() override { return _initial_guess; } + + T objective(const std::array& params) override { + return _objective_func(params); + } + + std::array gradient(const std::array& params) override { + return autoopt::gradient(_objective_func, params); + } + + matrix_t hessian(const std::array& params) override { + return autoopt::hessian(_objective_func, params); + } +}; + +template +struct log_barrier_optimization_problem + : public optimization_problem { + optimization_problem& _base_problem; + std::array _delta; + T _barrier_strength; + + log_barrier_optimization_problem( + optimization_problem& base_problem, + std::array delta, + T barrier_strength = T{1e-3}) + : _base_problem(base_problem), + _delta(delta), + _barrier_strength(barrier_strength) {} + + std::array initial_guess() override { + return _base_problem.initial_guess(); + } + + T objective(const std::array& params) override { + T base_obj = _base_problem.objective(params); + T barrier = barrier_term(params); + return base_obj + barrier; + } + + std::array gradient(const std::array& params) override { + auto base_grad = _base_problem.gradient(params); + std::array barrier_grad = autoopt::gradient( + [this](const std::array& p) { + return barrier_term(p); + }, + params); + std::array total_grad; + for (size_t i = 0; i < N; ++i) { + total_grad[i] = base_grad[i] + barrier_grad[i]; + } + return total_grad; + } + + matrix_t hessian(const std::array& params) override { + auto base_hess = _base_problem.hessian(params); + matrix_t barrier_hess = autoopt::hessian( + [this](const std::array& p) { + return barrier_term(p); + }, + params); + matrix_t total_hess; + for (size_t i = 0; i < N; ++i) { + for (size_t j = 0; j < N; ++j) { + total_hess[i][j] = base_hess[i][j] + barrier_hess[i][j]; + } + } + return total_hess; + } + +private: + template + U barrier_term(const std::array& params) { + U barrier = U{0}; + for (size_t i = 0; i < N; ++i) { + U lb = _base_problem.initial_guess()[i] - _delta[i]; + U ub = _base_problem.initial_guess()[i] + _delta[i]; + barrier = barrier + std::log(params[i] - lb) + std::log(ub - params[i]); + } + return -U{_barrier_strength} * barrier; + } +}; + +} // namespace autoopt diff --git a/include/autoopt/quadric.hpp b/include/autoopt/quadric.hpp new file mode 100644 index 0000000..3d8df50 --- /dev/null +++ b/include/autoopt/quadric.hpp @@ -0,0 +1,60 @@ +#pragma once + +#include +#include + +namespace autoopt { + +template +struct quadric { + T _A; // x² coefficient + T _B; // xy coefficient + T _C; // y² coefficient + T _D; // x coefficient + T _E; // y coefficient + T _F; // constant term + + constexpr quadric(T A = 0, T B = 0, T C = 0, T D = 0, T E = 0, T F = 0) + : _A(A), _B(B), _C(C), _D(D), _E(E), _F(F) {} + + constexpr quadric rotated_by(T angle_rad) const { + T co = std::cos(angle_rad); + T si = std::sin(angle_rad); + T co_2 = co * co; + T si_2 = si * si; + T co_si = co * si; + + T A_new = _A * co_2 - _B * co_si + _C * si_2; + T B_new = T{2} * (_A - _C) * co_si + _B * (co_2 - si_2); + T C_new = _A * si_2 + _B * co_si + _C * co_2; + T D_new = _D * co - _E * si; + T E_new = _D * si + _E * co; + T F_new = _F; + + return quadric(A_new, B_new, C_new, D_new, E_new, F_new); + } + + constexpr T at(T x) const { + T sign = T{-1}; + T Bx_E = _B * x + _E; + T Ax2_Dx_F = _A * x * x + _D * x + _F; + + if (_C == T{0}) { + return -Ax2_Dx_F / Bx_E; + } + + T discriminant = Bx_E * Bx_E - T{4} * _C * Ax2_Dx_F; + T num = sign * std::sqrt(discriminant) - Bx_E; + T denom = T{2} * _C; + return num / denom; + } + + constexpr T slope_at(T x) const { + return derivative([&](U x_val) { + quadric q{U{_A}, U{_B}, U{_C}, U{_D}, U{_E}, U{_F}}; + return q.at(x_val); + }, x); + } +}; + +} // namespace autoopt diff --git a/include/autoopt/util.hpp b/include/autoopt/util.hpp new file mode 100644 index 0000000..2853c97 --- /dev/null +++ b/include/autoopt/util.hpp @@ -0,0 +1,38 @@ +#pragma once + +#include + +namespace autoopt { + +template +T rad2deg(T radians) { + return radians * (T{180} / M_PI); +} + +template +T deg2rad(T degrees) { + return degrees * (M_PI / T{180}); +} + +template +T rad2arcmin(T radians) { + return radians * (T{180 * 60} / M_PI); +} + +template +T arcmin2rad(T arcminutes) { + return arcminutes * (M_PI / T{180 * 60}); +} + +template +T rad2arcsec(T radians) { + return radians * (T{180 * 60 * 60} / M_PI); +} + +template +T arcsec2rad(T arcseconds) { + return arcseconds * (M_PI / T{180 * 60 * 60}); +} + + +} // namespace autoopt diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 0000000..bbba31b --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,10 @@ +file(GLOB_RECURSE TEST_SOURCES *.cpp) + +add_executable(autoopt-test ${TEST_SOURCES}) +target_link_libraries(autoopt-test autoopt gtest gtest_main) + +install(TARGETS autoopt-test DESTINATION bin) + +include(GoogleTest) + +gtest_discover_tests(autoopt-test) \ No newline at end of file diff --git a/tests/dual.cpp b/tests/dual.cpp new file mode 100644 index 0000000..91cd555 --- /dev/null +++ b/tests/dual.cpp @@ -0,0 +1,137 @@ +#include "autoopt/dual.hpp" + +#include + +#include "autoopt/derivative.hpp" + +using namespace autoopt; + +TEST(DualTest, BasicOperations) { + dual a(2.0, 1.0); // a = 2.0, da/dx = 1.0 + dual b(3.0, 0.0); // b = 3.0, db/dx = 0.0 + + dual c = a + b; + EXPECT_DOUBLE_EQ(c._x, 5.0); + EXPECT_DOUBLE_EQ(c._dx, 1.0); + + dual d = a * b; + EXPECT_DOUBLE_EQ(d._x, 6.0); + EXPECT_DOUBLE_EQ(d._dx, 3.0); + + dual e = a / b; + EXPECT_DOUBLE_EQ(e._x, 2.0 / 3.0); + EXPECT_DOUBLE_EQ(e._dx, 1.0 / 3.0); +} + +TEST(DualTest, StandardFunctions) { + dual a(0.5, 1.0); // a = 0.5, da/dx = 1.0 + + dual b = std::sin(a); + EXPECT_DOUBLE_EQ(b._x, std::sin(0.5)); + EXPECT_DOUBLE_EQ(b._dx, std::cos(0.5)); + + dual c = std::exp(a); + EXPECT_DOUBLE_EQ(c._x, std::exp(0.5)); + EXPECT_DOUBLE_EQ(c._dx, std::exp(0.5)); + + dual d = std::log(a); + EXPECT_DOUBLE_EQ(d._x, std::log(0.5)); + EXPECT_DOUBLE_EQ(d._dx, 1.0 / 0.5); +} + +TEST(DualTest, DerivativeFunction) { + auto func = [](const T& x) { return std::sin(x) * std::exp(x); }; + + for (double val : {0.0, 0.5, 1.0, 2.0}) { + double deriv = derivative(func, val); + double expected = (std::cos(val) + std::sin(val)) * std::exp(val); + EXPECT_DOUBLE_EQ(deriv, expected); + } +} + +TEST(DualTest, GradientFunction) { + auto func = [](const std::array& x) { + return x[0] * x[0] + std::sin(x[1]); + }; + + std::array point = {1.0, 0.0}; + std::array grad = gradient(func, point); + + EXPECT_DOUBLE_EQ(grad[0], 2.0 * point[0]); // d/dx1 + EXPECT_DOUBLE_EQ(grad[1], std::cos(point[1])); // d/dx2 +} + +TEST(DualTest, JacobianFunction) { + auto func = [](const std::array& x) { + return std::array{x[0] * x[0], std::sin(x[1])}; + }; + + std::array point = {1.0, 0.0}; + auto jacob = jacobian(func, point); + + EXPECT_DOUBLE_EQ(jacob[0][0], 2.0 * point[0]); // d(f1)/d(x1) + EXPECT_DOUBLE_EQ(jacob[0][1], 0.0); // d(f1)/d(x2) + EXPECT_DOUBLE_EQ(jacob[1][0], 0.0); // d(f2)/d(x1) + EXPECT_DOUBLE_EQ(jacob[1][1], std::cos(point[1])); // d(f2)/d(x2) +} + +TEST(DualTest, HessianFunction) { + auto func = [](const std::array& x) { + return x[0] * x[0] + x[1] * x[1]; + }; + + std::array point = {1.0, 2.0}; + auto hess = hessian(func, point); + + EXPECT_DOUBLE_EQ(hess[0][0], 2.0); // d²f/dx1² + EXPECT_DOUBLE_EQ(hess[0][1], 0.0); // d²f/dx1dx2 + EXPECT_DOUBLE_EQ(hess[1][0], 0.0); // d²f/dx2dx1 + EXPECT_DOUBLE_EQ(hess[1][1], 2.0); // d²f/dx2² +} + +struct opti_func { + std::vector> test_data; + + // loss function + template + T operator()(const std::array& params) const { + T sum = T(0); + for (const auto& data_point : test_data) { + T x = T{data_point[0]}; + T y_true = T{data_point[1]}; + T a = params[0]; + T b = params[1]; + T c = params[2]; + T y_pred = a * x * x + b * x + c; + T error = y_pred - y_true; + sum = sum + error * error; + } + return sum / T(test_data.size()); + } +}; + +TEST(DualTest, OptimizationFunction) { + opti_func f; + f.test_data = { + {0.0, 4.0}, + {1.0, 1.0}, + {2.0, 0.0}, + {3.0, 1.0}, + {4.0, 4.0}, + }; + + std::array params = {1.0, -4.0, 4.0}; + + auto grad = gradient(f, params); + + EXPECT_DOUBLE_EQ(grad[0], 0.0); // dL/da + EXPECT_DOUBLE_EQ(grad[1], 0.0); // dL/db + EXPECT_DOUBLE_EQ(grad[2], 0.0); // dL/dc + + auto hess = hessian(f, params); + for (std::size_t i = 0; i < 3; ++i) { + for (std::size_t j = 0; j < 3; ++j) { + EXPECT_GE(hess[i][j], 0.0); // Hessian should be positive semi-definite + } + } +} diff --git a/tests/ellipse.cpp b/tests/ellipse.cpp new file mode 100644 index 0000000..41cd5dc --- /dev/null +++ b/tests/ellipse.cpp @@ -0,0 +1,78 @@ +#include + +#include +#include +#include +#include +#include + +using namespace autoopt; + +TEST(Ellipse, Slope) { + ellipse e{100, 1000, deg2rad(1.0)}; // entrance angle 1 degree + quadric q = e.to_quadric(); + + EXPECT_NEAR(q.slope_at(-10), -0.0010305116165301856, 1e-9); + EXPECT_NEAR(q.slope_at(0), 0.0, 1e-9); + EXPECT_NEAR(q.slope_at(10), 0.00090001261192696272, 1e-9); +} + +TEST(Ellipse, ParamGradient) { + std::vector> data_points = { + {-10.0, -0.001}, {0.0, 0.0}, {10, 0.0009}}; + + std::array params = {100, 1000, deg2rad(1.0), 0.0}; + + auto loss_func = [&data_points](const std::array& p) { + ellipse e{T{p[0]}, T{p[1]}, T{p[2]}}; + quadric q = e.to_quadric().rotated_by(T{p[3]}); + T loss = T{0}; + for (const auto& [x, y_true] : data_points) { + T y_pred = q.slope_at(T{x}); + T error = y_pred - T{y_true}; + loss = loss + error * error; + } + return loss / T(data_points.size()); + }; + + auto_diff_optimization_problem problem(loss_func, params); + + auto grad = problem.gradient(params); + EXPECT_NEAR(grad[0], -2.0789313126683308e-10, 1e-15); // d/d(left_arm) + EXPECT_NEAR(grad[1], -1.7464984353858657e-12, 1e-15); // d/d(right_arm) + EXPECT_NEAR(grad[2], 1.2013025455499119e-06, 1e-15); // d/d(entrance_angle) + EXPECT_NEAR(grad[3], -2.0332702665822054e-05, 1e-15); // d/d(rotation_angle) + + std::cout << "Gradient:\n"; + for (size_t i = 0; i < 4; ++i) { + std::cout << grad[i] << "\n"; + } + + auto hess = problem.hessian(params); + + // set formatting for easier reading + std::cout << std::scientific; + // set field width for alignment + std::cout << "Hessian matrix:\n"; + + for (size_t i = 0; i < 4; ++i) { + ; + for (size_t j = 0; j < 4; ++j) { + std::cout << std::setprecision(5) << std::setw(15) << hess[i][j]; + } + std::cout << "\n"; + } + + // log barrier + log_barrier_optimization_problem log_barrier_problem( + problem, + {1.0, 1.0, deg2rad(0.1), deg2rad(0.1)}, + 1e-3); + + auto log_barrier_grad = log_barrier_problem.gradient(params); + + std::cout << "Log Barrier Gradient:\n"; + for (size_t i = 0; i < 4; ++i) { + std::cout << log_barrier_grad[i] << "\n"; + } +} \ No newline at end of file diff --git a/tests/quadric.cpp b/tests/quadric.cpp new file mode 100644 index 0000000..bb58cb8 --- /dev/null +++ b/tests/quadric.cpp @@ -0,0 +1,21 @@ +#include + +#include "autoopt/quadric.hpp" +#include + +using namespace autoopt; + +TEST(QuadricTest, ParabolaRotation) { + quadric q(1.0, 0.0, 0.0, 0.0, -1.0, 0.0); // y = x^2 + + EXPECT_DOUBLE_EQ(q.at(-1.0), 1.0); // At x=1, y=1 + EXPECT_DOUBLE_EQ(q.at(0.0), 0.0); // At x=0, y=0 + EXPECT_DOUBLE_EQ(q.at(1.0), 1.0); // At x=1, y=1 + + double angle = M_PI / 4; // 45 degrees + quadric q_rotated = q.rotated_by(angle); + + EXPECT_NEAR(q_rotated.at(-1.0), -0.11729096183611623, 1e-9); + EXPECT_NEAR(q_rotated.at(0.0), 0.0, 1e-9); + EXPECT_TRUE(std::isnan(q_rotated.at(1.0))); // Expect NaN due to no real solution +} \ No newline at end of file