init

.gitignore

@@ -0,0 +1,9 @@
+# nix
+result
+
+# cmake
+build
+
+# direnv files
+.direnv
+.envrc

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()
+

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
+}

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"
+            '';
+          };
+        };
+      }
+    );
+}

include/autoopt/derivative.hpp

@@ -0,0 +1,61 @@
+#pragma once
+
+#include <array>
+
+#include "autoopt/dual.hpp"
+
+namespace autoopt {
+
+template <typename T, class Func>
+T derivative(Func&& f, const T& x) {
+  dual<T> a(x, T(1));
+  dual<T> b = f(a);
+  return b._dx;
+}
+
+template <typename T, std::size_t N, class Func>
+std::array<T, N> gradient(Func&& f, const std::array<T, N>& x) {
+  std::array<T, N> grad{};
+  std::array<dual<T>, N> dual_x{};
+  for (std::size_t i = 0; i < N; ++i) {
+    dual_x[i] = dual<T>(x[i], T(0));
+  }
+  for (std::size_t i = 0; i < N; ++i) {
+    dual_x[i]._dx = T(1);
+    dual<T> dual_y = f(dual_x);
+    grad[i] = dual_y._dx;
+    dual_x[i]._dx = T(0);
+  }
+  return grad;
+}
+
+template <typename T, std::size_t N, std::size_t M>
+using matrix_t = std::array<std::array<T, M>, N>;
+
+template <typename T, std::size_t N, std::size_t M, class Func>
+matrix_t<T, M, N> jacobian(Func&& f, const std::array<T, N>& x) {
+  matrix_t<T, M, N> jacob{};
+  std::array<dual<T>, N> dual_x{};
+  for (std::size_t i = 0; i < N; ++i) {
+    dual_x[i] = dual<T>(x[i], T(0));
+  }
+  for (std::size_t i = 0; i < N; ++i) {
+    dual_x[i]._dx = T(1);
+    std::array<dual<T>, 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 <typename T, std::size_t N, class Func>
+matrix_t<T, N, N> hessian(Func&& f, const std::array<T, N>& x) {
+  auto helper_func = [&f]<typename U>(const std::array<U, N>& y) {
+    return gradient<U, N>(f, y);
+  };
+  return jacobian<T, N, N>(helper_func, x);
+}
+
+} // namespace autoopt

include/autoopt/dual.hpp

@@ -0,0 +1,241 @@
+#pragma once
+
+#include <cmath>
+
+namespace autoopt {
+
+template <typename T>
+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 <typename U>
+  constexpr dual(U x) : _x(T(x)), _dx(T(0)) {}
+};
+
+template <typename T>
+constexpr dual<T> operator+(const dual<T>& a, const dual<T>& b) {
+  return dual<T>(a._x + b._x, a._dx + b._dx);
+}
+
+template <typename T>
+constexpr dual<T> operator-(const dual<T>& a) {
+  return dual<T>(-a._x, -a._dx);
+}
+
+template <typename T>
+constexpr dual<T> operator-(const dual<T>& a, const dual<T>& b) {
+  return dual<T>(a._x - b._x, a._dx - b._dx);
+}
+
+template <typename T>
+constexpr dual<T> operator*(const dual<T>& a, const dual<T>& b) {
+  return dual<T>(a._x * b._x, a._x * b._dx + a._dx * b._x);
+}
+
+template <typename T>
+constexpr dual<T> operator/(const dual<T>& a, const dual<T>& b) {
+  return dual<T>(a._x / b._x, (a._dx * b._x - a._x * b._dx) / (b._x * b._x));
+}
+
+template <typename T>
+constexpr bool operator==(const dual<T>& a, const dual<T>& b) {
+  return a._x == b._x;
+}
+
+template <typename T>
+constexpr bool operator!=(const dual<T>& a, const dual<T>& b) {
+  return a._x != b._x;
+}
+
+template <typename T>
+constexpr bool operator<(const dual<T>& a, const dual<T>& b) {
+  return a._x < b._x;
+}
+
+template <typename T>
+constexpr bool operator<=(const dual<T>& a, const dual<T>& b) {
+  return a._x <= b._x;
+}
+
+template <typename T>
+constexpr bool operator>(const dual<T>& a, const dual<T>& b) {
+  return a._x > b._x;
+}
+
+template <typename T>
+constexpr bool operator>=(const dual<T>& a, const dual<T>& b) {
+  return a._x >= b._x;
+}
+
+}  // namespace autoopt
+
+namespace std {
+
+using autoopt::dual;
+
+// forward declarations of standard functions
+
+template <typename T>
+constexpr dual<T> abs(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> exp(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> log(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> pow(const dual<T>& a, const dual<T>& b);
+
+template <typename T>
+constexpr dual<T> sqrt(const dual<T>& a);
+
+// forward declarations of trigonometric functions
+
+template <typename T>
+constexpr dual<T> sin(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> cos(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> tan(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> asin(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> acos(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> atan(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> atan2(const dual<T>& y, const dual<T>& x);
+
+// forward declarations of hyperbolic functions
+
+template <typename T>
+constexpr dual<T> sinh(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> cosh(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> tanh(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> asinh(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> acosh(const dual<T>& a);
+
+template <typename T>
+constexpr dual<T> atanh(const dual<T>& a);
+
+// standard functions
+
+template <typename T>
+constexpr dual<T> abs(const dual<T>& a) {
+  return dual<T>(std::abs(a._x), (a._x >= T(0) ? T(1) : T(-1)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> exp(const dual<T>& a) {
+  T exp_x = std::exp(a._x);
+  return dual<T>(exp_x, exp_x * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> log(const dual<T>& a) {
+  return dual<T>(std::log(a._x), (T(1) / a._x) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> pow(const dual<T>& a, const dual<T>& b) {
+  return std::exp(b * std::log(a));
+}
+
+template <typename T>
+constexpr dual<T> sqrt(const dual<T>& a) {
+  T sqrt_x = std::sqrt(a._x);
+  return dual<T>(sqrt_x, (T(1) / (T(2) * sqrt_x)) * a._dx);
+}
+
+// trigonometric functions
+
+template <typename T>
+constexpr dual<T> sin(const dual<T>& a) {
+  return dual<T>(std::sin(a._x), std::cos(a._x) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> cos(const dual<T>& a) {
+  return dual<T>(std::cos(a._x), -std::sin(a._x) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> tan(const dual<T>& a) {
+  T cos_x = std::cos(a._x);
+  return dual<T>(std::tan(a._x), (T(1) / (cos_x * cos_x)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> asin(const dual<T>& a) {
+  return dual<T>(std::asin(a._x), (T(1) / std::sqrt(T(1) - a._x * a._x)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> acos(const dual<T>& a) {
+  return dual<T>(std::acos(a._x), (-T(1) / std::sqrt(T(1) - a._x * a._x)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> atan(const dual<T>& a) {
+  return dual<T>(std::atan(a._x), (T(1) / (T(1) + a._x * a._x)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> atan2(const dual<T>& y, const dual<T>& x) {
+  return dual<T>(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 <typename T>
+constexpr dual<T> sinh(const dual<T>& a) {
+  return dual<T>(std::sinh(a._x), std::cosh(a._x) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> cosh(const dual<T>& a) {
+  return dual<T>(std::cosh(a._x), std::sinh(a._x) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> tanh(const dual<T>& a) {
+  T cosh_x = std::cosh(a._x);
+  return dual<T>(std::tanh(a._x), (T(1) / (cosh_x * cosh_x)) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> asinh(const dual<T>& a) {
+  return dual<T>(std::asinh(a._x), (T(1) / std::sqrt(a._x * a._x + T(1))) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> acosh(const dual<T>& a) {
+  return dual<T>(std::acosh(a._x), (T(1) / std::sqrt(a._x * a._x - T(1))) * a._dx);
+}
+
+template <typename T>
+constexpr dual<T> atanh(const dual<T>& a) {
+  return dual<T>(std::atanh(a._x), (T(1) / (T(1) - a._x * a._x)) * a._dx);
+}
+
+}  // namespace std

include/autoopt/ellipse.hpp

@@ -0,0 +1,50 @@
+#pragma once
+
+#include "autoopt/quadric.hpp"
+
+namespace autoopt {
+
+template <typename T>
+struct ellipse {
+  T left_arm;
+  T right_arm;
+  T entrance_angle;
+
+  constexpr quadric<T> 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<T>(A, B, C, D, E, F);
+  }
+};
+
+} // namespace autoopt

include/autoopt/optimization_problem.hpp

@@ -0,0 +1,111 @@
+#pragma once
+
+#include <array>
+#include <cstddef>
+
+#include "autoopt/derivative.hpp"
+
+namespace autoopt {
+
+template <typename T, size_t N>
+struct optimization_problem {
+  virtual std::array<T, N> initial_guess() = 0;
+
+  virtual T objective(const std::array<T, N>& params) = 0;
+  virtual std::array<T, N> gradient(const std::array<T, N>& params) = 0;
+  virtual matrix_t<T, N, N> hessian(const std::array<T, N>& params) = 0;
+};
+
+template <typename T, size_t N, class Func>
+struct auto_diff_optimization_problem : public optimization_problem<T, N> {
+  Func _objective_func;
+  std::array<T, N> _initial_guess;
+
+  auto_diff_optimization_problem(
+      Func objective_func, std::array<T, N> initial_guess = std::array<T, N>{})
+      : _objective_func(objective_func), _initial_guess(initial_guess) {}
+
+  std::array<T, N> initial_guess() override { return _initial_guess; }
+
+  T objective(const std::array<T, N>& params) override {
+    return _objective_func(params);
+  }
+
+  std::array<T, N> gradient(const std::array<T, N>& params) override {
+    return autoopt::gradient<T, N>(_objective_func, params);
+  }
+
+  matrix_t<T, N, N> hessian(const std::array<T, N>& params) override {
+    return autoopt::hessian<T, N>(_objective_func, params);
+  }
+};
+
+template <typename T, size_t N>
+struct log_barrier_optimization_problem
+    : public optimization_problem<T, N> {
+  optimization_problem<T, N>& _base_problem;
+  std::array<T, N> _delta;
+  T _barrier_strength;
+
+  log_barrier_optimization_problem(
+      optimization_problem<T, N>& base_problem,
+      std::array<T, N> delta,
+      T barrier_strength = T{1e-3})
+      : _base_problem(base_problem),
+        _delta(delta),
+        _barrier_strength(barrier_strength) {}
+
+  std::array<T, N> initial_guess() override {
+    return _base_problem.initial_guess();
+  }
+
+  T objective(const std::array<T, N>& params) override {
+    T base_obj = _base_problem.objective(params);
+    T barrier = barrier_term(params);
+    return base_obj + barrier;
+  }
+
+  std::array<T, N> gradient(const std::array<T, N>& params) override {
+    auto base_grad = _base_problem.gradient(params);
+    std::array<T, N> barrier_grad = autoopt::gradient<T, N>(
+        [this]<typename U>(const std::array<U, N>& p) {
+          return barrier_term<U>(p);
+        },
+        params);
+    std::array<T, N> total_grad;
+    for (size_t i = 0; i < N; ++i) {
+      total_grad[i] = base_grad[i] + barrier_grad[i];
+    }
+    return total_grad;
+  }
+
+  matrix_t<T, N, N> hessian(const std::array<T, N>& params) override {
+    auto base_hess = _base_problem.hessian(params);
+    matrix_t<T, N, N> barrier_hess = autoopt::hessian<T, N>(
+        [this]<typename U>(const std::array<U, N>& p) {
+          return barrier_term<U>(p);
+        },
+        params);
+    matrix_t<T, N, N> 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 <typename U>
+  U barrier_term(const std::array<U, N>& 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

include/autoopt/quadric.hpp

@@ -0,0 +1,60 @@
+#pragma once
+
+#include <cmath>
+#include <autoopt/derivative.hpp>
+
+namespace autoopt {
+
+template <typename T>
+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([&]<typename U>(U x_val) {
+      quadric<U> q{U{_A}, U{_B}, U{_C}, U{_D}, U{_E}, U{_F}};
+      return q.at(x_val);
+    }, x);
+  }
+};
+
+}  // namespace autoopt

include/autoopt/util.hpp

@@ -0,0 +1,38 @@
+#pragma once
+
+#include <cmath>
+
+namespace autoopt {
+
+template <typename T>
+T rad2deg(T radians) {
+  return radians * (T{180} / M_PI);
+}
+
+template <typename T>
+T deg2rad(T degrees) {
+  return degrees * (M_PI / T{180}); 
+}
+
+template <typename T>
+T rad2arcmin(T radians) {
+  return radians * (T{180 * 60} / M_PI);
+}
+
+template <typename T>
+T arcmin2rad(T arcminutes) {
+  return arcminutes * (M_PI / T{180 * 60});
+}
+
+template <typename T>
+T rad2arcsec(T radians) {
+  return radians * (T{180 * 60 * 60} / M_PI);
+}
+
+template <typename T>
+T arcsec2rad(T arcseconds) {
+  return arcseconds * (M_PI / T{180 * 60 * 60});
+}
+
+
+}  // namespace autoopt

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)

tests/dual.cpp

@@ -0,0 +1,137 @@
+#include "autoopt/dual.hpp"
+
+#include <gtest/gtest.h>
+
+#include "autoopt/derivative.hpp"
+
+using namespace autoopt;
+
+TEST(DualTest, BasicOperations) {
+  dual<double> a(2.0, 1.0);  // a = 2.0, da/dx = 1.0
+  dual<double> b(3.0, 0.0);  // b = 3.0, db/dx = 0.0
+
+  dual<double> c = a + b;
+  EXPECT_DOUBLE_EQ(c._x, 5.0);
+  EXPECT_DOUBLE_EQ(c._dx, 1.0);
+
+  dual<double> d = a * b;
+  EXPECT_DOUBLE_EQ(d._x, 6.0);
+  EXPECT_DOUBLE_EQ(d._dx, 3.0);
+
+  dual<double> 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<double> a(0.5, 1.0);  // a = 0.5, da/dx = 1.0
+
+  dual<double> b = std::sin(a);
+  EXPECT_DOUBLE_EQ(b._x, std::sin(0.5));
+  EXPECT_DOUBLE_EQ(b._dx, std::cos(0.5));
+
+  dual<double> c = std::exp(a);
+  EXPECT_DOUBLE_EQ(c._x, std::exp(0.5));
+  EXPECT_DOUBLE_EQ(c._dx, std::exp(0.5));
+
+  dual<double> 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 = []<typename T>(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 = []<typename T>(const std::array<T, 2>& x) {
+    return x[0] * x[0] + std::sin(x[1]);
+  };
+
+  std::array<double, 2> point = {1.0, 0.0};
+  std::array<double, 2> 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 = []<typename T>(const std::array<T, 2>& x) {
+    return std::array<T, 2>{x[0] * x[0], std::sin(x[1])};
+  };
+
+  std::array<double, 2> point = {1.0, 0.0};
+  auto jacob = jacobian<double, 2, 2>(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 = []<typename T>(const std::array<T, 2>& x) {
+    return x[0] * x[0] + x[1] * x[1];
+  };
+
+  std::array<double, 2> point = {1.0, 2.0};
+  auto hess = hessian<double, 2>(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<std::array<double, 2>> test_data;
+
+  // loss function
+  template <typename T>
+  T operator()(const std::array<T, 3>& 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<double, 3> params = {1.0, -4.0, 4.0};
+  
+  auto grad = gradient<double, 3>(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<double, 3>(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
+    }
+  }
+}

tests/ellipse.cpp

@@ -0,0 +1,78 @@
+#include <gtest/gtest.h>
+
+#include <autoopt/ellipse.hpp>
+#include <autoopt/optimization_problem.hpp>
+#include <autoopt/util.hpp>
+#include <iomanip>
+#include <iostream>
+
+using namespace autoopt;
+
+TEST(Ellipse, Slope) {
+  ellipse<double> e{100, 1000, deg2rad(1.0)};  // entrance angle 1 degree
+  quadric<double> 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<std::pair<double, double>> data_points = {
+      {-10.0, -0.001}, {0.0, 0.0}, {10, 0.0009}};
+
+  std::array<double, 4> params = {100, 1000, deg2rad(1.0), 0.0};
+
+  auto loss_func = [&data_points]<typename T>(const std::array<T, 4>& p) {
+    ellipse<T> e{T{p[0]}, T{p[1]}, T{p[2]}};
+    quadric<T> 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<double, 4, decltype(loss_func)> 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<double, 4> 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";
+  }
+}

tests/quadric.cpp

@@ -0,0 +1,21 @@
+#include <gtest/gtest.h>
+
+#include "autoopt/quadric.hpp"
+#include <iostream>
+
+using namespace autoopt;
+
+TEST(QuadricTest, ParabolaRotation) {
+  quadric<double> 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<double> 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
+}