add blts

include/autoopt/derivative.hpp

@@ -1,6 +1,6 @@
 #pragma once
 
-#include <array>
+#include <Eigen/Eigen>
 
 #include "autoopt/dual.hpp"
 
@@ -13,49 +13,46 @@ T derivative(Func&& f, const T& x) {
   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));
+template <typename T, class Func>
+Eigen::VectorX<T> gradient(Func&& f, const Eigen::VectorX<T>& x) {
+  Eigen::VectorX<T> grad{x.size()};
+  Eigen::VectorX<dual<T>> dual_x{x.size()};
+  for (int i = 0; i < x.size(); ++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);
+  for (int i = 0; i < x.size(); ++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);
+    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));
+template <typename T, class Func>
+Eigen::MatrixX<T> jacobian(Func&& f, const Eigen::VectorX<T>& x) {
+  Eigen::MatrixX<T> jacob(f(x).size(), x.size());
+  Eigen::VectorX<dual<T>> dual_x(x.size());
+  for (int i = 0; i < x.size(); ++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;
+  for (int i = 0; i < x.size(); ++i) {
+    dual_x(i)._dx = T(1);
+    Eigen::VectorX<dual<T>> dual_y = f(dual_x);
+    for (int j = 0; j < dual_y.size(); ++j) {
+      jacob(j, i) = dual_y(j)._dx;
     }
-    dual_x[i]._dx = T(0);
+    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);
+template <typename T, class Func>
+Eigen::MatrixX<T> hessian(Func&& f, const Eigen::VectorX<T>& x) {
+  auto helper_func = [&f]<typename U>(const Eigen::VectorX<U>& y) {
+    return gradient<U>(f, y);
   };
-  return jacobian<T, N, N>(helper_func, x);
+  return jacobian<T>(helper_func, x);
 }
 
 } // namespace autoopt