rework

.clang-format

@@ -0,0 +1,3 @@
+---
+BasedOnStyle: Google
+IndentWidth: 2

.gitignore

@@ -4,6 +4,9 @@ result
 # cmake
 build
 
+# clangd
+.cache
+
 # direnv files
 .direnv
 .envrc

CMakeLists.txt

@@ -1,23 +1,48 @@
 cmake_minimum_required(VERSION 4.0)
 
-project(autoopt)
-
+set(CMAKE_C_COMPILER clang)
+set(CMAKE_CXX_COMPILER clang++)
 set(CMAKE_CXX_STANDARD 23)
 
+project(slopefit)
+
 add_compile_options(-Wall -Werror -Wpedantic)
 
 find_package(Eigen3 REQUIRED)
 
-add_library(autoopt INTERFACE)
-target_include_directories(autoopt INTERFACE include)
-target_link_libraries(autoopt INTERFACE Eigen3::Eigen)
+find_package(Python 3.13 COMPONENTS Interpreter Development.Module REQUIRED)
 
-install(DIRECTORY include/autoopt DESTINATION include)
+execute_process(
+  COMMAND "${Python_EXECUTABLE}" -m nanobind --cmake_dir
+  OUTPUT_STRIP_TRAILING_WHITESPACE OUTPUT_VARIABLE nanobind_ROOT)
+find_package(nanobind CONFIG REQUIRED)
+
+file(GLOB_RECURSE SOURCES lib/*.cpp)
+
+add_library(autoopt STATIC ${SOURCES})
+target_include_directories(autoopt PUBLIC include)
+target_link_libraries(autoopt PUBLIC Eigen3::Eigen)
+target_link_libraries(autoopt PRIVATE fftw3)
+
+nanobind_add_module(slopefit NOMINSIZE src/module.cpp)
+nanobind_add_stub(
+    slopefit_stub
+    MODULE slopefit
+    OUTPUT slopefit.pyi
+    PYTHON_PATH $<TARGET_FILE_DIR:slopefit>
+    DEPENDS slopefit
+)
+
+target_link_libraries(slopefit PRIVATE autoopt)
+install(TARGETS slopefit DESTINATION .)
+install(FILES $<TARGET_FILE_DIR:slopefit>/slopefit.pyi DESTINATION .)
 
 # add tests if gtest is found
+enable_testing()
+add_subdirectory(tests)
+
 find_library(GTestPackage gtest QUIET)
 if(GTestPackage)
-  enable_testing()
-  add_subdirectory(tests)
+
 endif()
 

README.md

@@ -0,0 +1 @@
+# slopefit

build.nix

@@ -0,0 +1,44 @@
+{
+  cmake,
+  ninja,
+  eigen,
+  fftw,
+  clang-tools,
+  python3Packages,
+  doCheck ? false,
+  gtest,
+}:
+python3Packages.buildPythonPackage {
+  pname = "slopefit";
+  version = "0.1.0";
+  pyproject = true;
+
+  dependencies = [
+    python3Packages.numpy
+    python3Packages.matplotlib
+    python3Packages.pyside6
+    python3Packages.pandas
+    python3Packages.pyqtgraph
+  ];
+
+  buildInputs = [
+    eigen
+    fftw
+  ]
+  ++ (if doCheck then [ gtest ] else [ ]);
+
+  nativeBuildInputs = [
+    cmake
+    clang-tools
+    ninja
+    python3Packages.scikit-build-core
+    python3Packages.nanobind
+  ];
+
+  buildPhase = ''
+    cd ..
+    pypaBuildPhase
+  '';
+
+  src = ./.;
+}

flake.lock

@@ -20,11 +20,11 @@
     },
     "nixpkgs": {
       "locked": {
-        "lastModified": 1768564909,
-        "narHash": "sha256-Kell/SpJYVkHWMvnhqJz/8DqQg2b6PguxVWOuadbHCc=",
+        "lastModified": 1789546076,
+        "narHash": "sha256-zVxLZiSnmaaPLwnhj7pwmqe3axBg/C6nG5JZsJMh2g4=",
         "owner": "nixos",
         "repo": "nixpkgs",
-        "rev": "e4bae1bd10c9c57b2cf517953ab70060a828ee6f",
+        "rev": "b1b875982b17dabde9b4a37f3e229e74913e6db3",
         "type": "github"
       },
       "original": {

flake.nix

@@ -1,5 +1,5 @@
 {
-  description = "auto-opt-cpp";
+  description = "slopefit";
 
   inputs = {
     nixpkgs.url = "github:nixos/nixpkgs?ref=nixos-unstable";
@@ -16,55 +16,24 @@
       system:
       let
         pkgs = nixpkgs.legacyPackages.${system};
-        pkg =
-          {
-            stdenv,
-            cmake,
-            ninja,
-            lib,
-            eigen,
-            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
-              eigen
-            ];
-          };
+        pkg = import ./build.nix;
       in
       rec {
         packages = rec {
-          autoopt-cpp = pkgs.callPackage pkg { };
-          autoopt-cpp-with-tests = pkgs.callPackage pkg { withTests = true; };
-          default = autoopt-cpp;
+          slopefit = pkgs.callPackage pkg { };
+          slopefit-with-tests = pkgs.callPackage pkg { doCheck = true; };
+          default = slopefit;
         };
         devShells = {
           default = pkgs.mkShell {
             inputsFrom = [
-              packages.autoopt-cpp-with-tests
+              packages.slopefit-with-tests
             ];
 
+            buildInputs = [ pkgs.python3Packages.pip ];
+
             shellHook = ''
-              echo "Development shell for autoopt-cpp"
+              echo "Development shell for slopefit"
             '';
           };
         };

include/autoopt/btls.hpp

@@ -1,7 +1,5 @@
 #pragma once
 
-#include <iostream>
-
 #include "autoopt/optimization_problem.hpp"
 
 namespace autoopt {
@@ -9,10 +7,10 @@ namespace autoopt {
 template <typename T>
 struct btls_parameters {
   T step_decrease = T{0.5};
-  T step_increase = T{1.5};
+  T step_increase = T{1.2};
   T sufficient_decrease = T{1e-2};
-  T tolerance = T{1e-9};
-  size_t max_iters = 1000;
+  T tolerance = T{1e-10};
+  size_t max_iters = 2000;
 };
 
 template <typename T>
@@ -24,20 +22,27 @@ void btls(optimization_problem<T>& problem,
   for (size_t iter = 0; iter < params.max_iters; ++iter) {
     T obj_value = problem.objective(x);
 
-    std::cout << "Iter " << iter << ": obj = " << obj_value
-              << ", x = " << x.transpose() << ", step_size = " << step_size
-              << std::endl;
+    Eigen::VectorX<T> grad = problem.gradient(x);
 
-    Eigen::VectorX<T> grad = -problem.gradient(x);
+    Eigen::MatrixX<T> hess = problem.hessian(x);
+    Eigen::VectorX<T> step_dir = -hess.ldlt().solve(grad).normalized();
 
-    Eigen::VectorX<T> step_dir = grad.normalized();
-    while (problem.objective(x + step_size * step_dir) >
-           obj_value +
-               params.sufficient_decrease * step_size * grad.dot(step_dir)) {
+    // Eigen::VectorX<T> step_dir = -grad.normalized();
+
+    auto decrease_condition = [&] {
+      T new_obj = problem.objective(x + step_size * step_dir);
+      return std::isnan(new_obj) ||
+             (new_obj > obj_value - std::abs(params.sufficient_decrease *
+                                             step_size * grad.dot(step_dir)));
+    };
+
+    while (decrease_condition()) {
       step_size *= params.step_decrease;
     }
     x += step_size * step_dir;
 
+    step_size = step_size * params.step_increase;
+
     if (step_size < params.tolerance) {
       break;
     }

include/autoopt/ellipse_new.hpp

@@ -1,16 +1,25 @@
 #pragma once
 
+#include "autoopt/parameters.hpp"
 #include "autoopt/quadric.hpp"
 
 namespace autoopt {
 
-template <typename T>
-struct ellipse {
-  T left_arm;
-  T right_arm;
-  T entrance_angle;
+struct ellipse_new : curve_base {
+  // T focal_length;
+  ellipse_new() : curve_base() {
+    initial = Eigen::VectorXd::Zero(4);
+    delta = Eigen::VectorXd::Zero(4);
+    param_names =
+        std::vector<const char*>{"left arm", "right arm", "entrance angle", "rotation"};
+  }
 
-  constexpr quadric<T> to_quadric() const {
+  template <typename T>
+  quadric<T> grapher_from_params(const Eigen::VectorX<T>& params) const {
+    T left_arm = params(0);
+    T right_arm = params(1);
+    T entrance_angle = params(2);
+    T alpha = params(3);
     T a = (left_arm + right_arm) / T{2};
     T a2 = a * a;
     
@@ -43,7 +52,7 @@ struct ellipse {
     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);
+    return quadric<T>(A, B, C, D, E, F).rotated_by(alpha);
   }
 };
 

include/autoopt/importer.hpp

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <vector>
+#include <Eigen/Eigen>
+#include <generator>
+#include <filesystem>
+#include <fstream>
+
+namespace autoopt {
+struct raw_data {
+  Eigen::VectorXd _xs;
+  std::vector<Eigen::VectorXd> _yss;
+};
+
+std::generator<std::pair<double, double>> parse_dat_file(
+    const std::filesystem::path& filename);
+
+std::pair<std::vector<double>, std::vector<double>> parse_row(auto& iter,
+                                                              auto sentinel);
+
+raw_data read_data(const std::filesystem::path& filename);
+
+}

include/autoopt/optimization_problem.hpp

@@ -1,7 +1,5 @@
 #pragma once
 
-#include <cstddef>
-
 #include "autoopt/derivative.hpp"
 
 namespace autoopt {
@@ -76,7 +74,7 @@ struct log_barrier_optimization_problem
     auto base_grad = _base_problem.gradient(params);
     Eigen::VectorX<T> barrier_grad = autoopt::gradient<T>(
         [this]<typename U>(const Eigen::VectorX<U>& p) {
-          return barrier_term<U>(p);
+          return this->barrier_term<U>(p);
         },
         params);
     Eigen::VectorX<T> total_grad(params.size());
@@ -88,7 +86,7 @@ struct log_barrier_optimization_problem
     auto base_hess = _base_problem.hessian(params);
     Eigen::MatrixX<T> barrier_hess = autoopt::hessian<T>(
         [this]<typename U>(const Eigen::VectorX<U>& p) {
-          return barrier_term<U>(p);
+          return this->barrier_term<U>(p);
         },
         params);
     Eigen::MatrixX<T> total_hess(params.size(), params.size());

include/autoopt/optimizer.hpp

@@ -0,0 +1,64 @@
+#pragma once
+
+#include <filesystem>
+
+#include <functional>
+#include <variant>
+#include <vector>
+#include "autoopt/parabola_new.hpp"
+#include "autoopt/ellipse_new.hpp"
+#include "autoopt/importer.hpp"
+
+namespace autoopt {
+
+enum CurveType {
+  Parabola,
+  Ellipse
+};
+
+using curves = std::variant<parabola_new, ellipse_new>;
+
+struct optimizer {
+
+public:
+  void import_file(std::filesystem::path file_path);
+  void set_lower_bound(double lb);
+  void set_upper_bound(double ub);
+
+  void set_curve_type(int curve_type);
+  void set_initial(size_t i, double value);
+  void set_delta(size_t i, double value);
+  void flip_curve();
+  void calc_fit();
+
+  void set_plot_input_cb(std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)> cb);
+  void set_plot_slope_cb(std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys, Eigen::VectorXd ys_curve)> cb);
+  void set_init_bounds_cb(std::function<void(double lb, double ub, double step)> cb);
+  void set_update_table_cb(std::function<void(std::vector<std::tuple<std::string, double, double>> entries)> cb);
+  void set_plot_error_cb(std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> cb);
+
+private:
+  struct data {
+    raw_data _raw;
+    double lb, ub;
+  };
+  std::optional<data> _data;
+  std::optional<test_data> _test_data;
+  curves _curve;
+  std::optional<std::pair<Eigen::VectorXd, double>> _fit_result;
+
+  std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)> _plot_input_cb;
+  std::function<void(double lb, double ub, double step)> _init_bounds_cb;
+  std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys, Eigen::VectorXd ys_curve)> _plot_slope_cb;
+  std::function<void(std::vector<std::tuple<std::string, double, double>> entries)> _param_table_cb;
+  std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> _plot_error_cb;
+
+  void update_test_data();
+  void init_bounds_cb();
+  void plot_input_cb();
+  void plot_slope_cb();
+  void update_table_cb();
+  void plot_error_cb();
+};
+
+}  // namespace autoopt

include/autoopt/parabola_new.hpp

@@ -0,0 +1,33 @@
+#pragma once
+
+#include "autoopt/parameters.hpp"
+#include "autoopt/quadric.hpp"
+
+namespace autoopt {
+
+struct parabola_new : curve_base {
+  // T focal_length;
+  parabola_new() : curve_base() {
+    initial = Eigen::VectorXd::Zero(3);
+    delta = Eigen::VectorXd::Zero(3);
+    param_names =
+        std::vector<const char*>{"focus", "entrance angle", "rotation"};
+  }
+
+  template <typename T>
+  quadric<T> grapher_from_params(const Eigen::VectorX<T>& params) const {
+    T exit_arm = params(0);
+    T entrance_angle = params(1);
+    T alpha = params(2);
+
+    T x = exit_arm * std::sin(T{2.0} * entrance_angle);
+    T f = T{0.5} * (T{1.0} - std::cos(T{2.0} * entrance_angle)) * exit_arm;
+    T a = T{1.0} / (T{4.0} * f);
+    T y = a * x * x;
+    quadric<T> q = quadric(a, T{0}, T{0}, T{0}, T{-1.0}, T{0});
+    return q.translated_by(-x, -y).rotated_by(entrance_angle - T{M_PI_2} +
+                                              alpha);
+  }
+};
+
+}  // namespace autoopt

include/autoopt/parameters.hpp

@@ -0,0 +1,76 @@
+#pragma once
+
+#include <Eigen/Eigen>
+#include <autoopt/btls.hpp>
+#include <autoopt/optimization_problem.hpp>
+#include <cmath>
+#include <cstddef>
+#include <optional>
+#include <utility>
+#include <vector>
+
+#include "autoopt/util.hpp"
+
+namespace autoopt {
+
+struct test_data {
+  Eigen::VectorXd _xs;
+  Eigen::VectorXd _ys;
+  double _center;
+
+  size_t size() const { return _xs.size(); }
+
+  std::pair<double, double> operator[](size_t i) const {
+    return {_xs(i) - _center, _ys(i)};
+  };
+};
+
+struct curve_base {
+  Eigen::VectorXd initial;
+  Eigen::VectorXd delta;
+  std::vector<const char*> param_names;
+  bool flip = false;
+
+  std::pair<Eigen::VectorXd, double> fit(this auto&& self,
+                                         const test_data& measurements) {
+    auto opt_func = [&]<typename T>(const Eigen::VectorX<T>& params) {
+      auto grapher = self.template grapher_from_params<T>(params);
+      T error{0.0};
+      for (size_t i = 0; i < measurements.size(); i++) {
+        auto [x, y] = measurements[i];
+        T distance = T{self.flip ? -1.0 : 1.0} *
+                         std::atan(grapher.slope_at(T{self.flip ? -x : x})) -
+                     T{arcsec2rad(y)};
+        error = error + distance * distance;
+      }
+      return error / T(measurements.size());
+    };
+
+    auto_diff_optimization_problem problem(opt_func, self.initial);
+
+    log_barrier_optimization_problem lb_problem(problem, self.delta, 1e-4);
+
+    while (lb_problem._barrier_strength > 1e-20) {
+      btls(lb_problem);
+      lb_problem._barrier_strength *= 1e-2;
+    }
+
+    Eigen::VectorX<double> fitted_params = problem.x();
+    return {fitted_params, std::sqrt(problem.objective(fitted_params))};
+  }
+
+  Eigen::VectorXd plot(this auto&& self, const test_data& measurements,
+                       std::optional<Eigen::VectorXd> params = {}) {
+    auto grapher = self.template grapher_from_params<double>(
+        params.value_or(self.initial));
+    Eigen::VectorXd ys(measurements.size());
+    for (size_t i = 0; i < measurements.size(); i++) {
+      double x = measurements[i].first;
+      ys(i) = rad2arcsec((self.flip ? -1.0 : 1.0) *
+                         atan(grapher.slope_at(self.flip ? -x : x)));
+    }
+    return ys;
+  }
+};
+
+}  // namespace autoopt

include/autoopt/quadric.hpp

@@ -1,7 +1,7 @@
 #pragma once
 
-#include <cmath>
 #include <autoopt/derivative.hpp>
+#include <cmath>
 
 namespace autoopt {
 
@@ -34,6 +34,14 @@ struct quadric {
     return quadric(A_new, B_new, C_new, D_new, E_new, F_new);
   }
 
+  constexpr quadric translated_by(T x_o, T y_o) const {
+    T D_new = _D - T{2} * _A * x_o - _B * y_o;
+    T E_new = _E - _B * x_o - T{2} * _C * y_o;
+    T F_new = _F + _A * x_o * x_o + _B * x_o * y_o + _C * y_o * y_o - _D * x_o -
+              _E * y_o;
+    return quadric(_A, _B, _C, D_new, E_new, F_new);
+  }
+
   constexpr T at(T x) const {
     T sign = T{-1};
     T Bx_E = _B * x + _E;
@@ -50,10 +58,12 @@ struct quadric {
   }
 
   constexpr T slope_at(T x) const {
-    return derivative([&]<typename U>(U x_val) {
+    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);
+        },
+        x);
   }
 };
 

include/autoopt/reactive.hpp

@@ -0,0 +1,79 @@
+
+#include <array>
+#include <cstdint>
+#include <functional>
+#include <tuple>
+
+namespace autoopt {
+
+template <typename T>
+struct source {
+  T _value;
+  std::uint64_t _version = 0;
+
+  source() = default;
+  explicit source(T value) : _value(std::move(value)) {}
+
+  virtual T& get() { return _value; };
+
+  bool set(T value) {
+    if (value == _value) return false;
+    _value = std::move(value);
+    ++_version;
+    return true;
+  }
+};
+
+template <typename O, typename... Is>
+struct transform : source<O> {
+  std::function<O(Is&...)> _func;
+  std::tuple<source<Is>&...> _sources;
+  std::array<std::uint64_t, sizeof...(Is)> _versions{};
+
+  transform(std::function<O(Is&...)> func, source<Is>&... sources)
+      : _func(func), _sources(sources...) {
+    this->_value = std::apply([&](source<Is>&... sources) {
+      return _func(sources.get()...);
+    }, _sources);
+    this->_versions = current_versions();
+  }
+
+  virtual O& get() override {
+    auto cur_versions = current_versions();
+    if (_versions == cur_versions) return this->_value;
+
+    O next = std::apply([&](source<Is>&... sources) {
+      return _func(sources.get()...);
+    }, _sources);
+
+    if (next != this->_value) {
+      this->_value = next;
+      ++this->_version;
+    }
+    return this->_value;
+  }
+
+ private: 
+  std::array<std::uint64_t, sizeof...(Is)> current_versions() {
+    return std::apply([&](source<Is>&... sources){
+      
+      return std::array<std::uint64_t, sizeof...(Is)>{sources._version...};
+    }, _sources);
+  }
+};
+
+template <typename T>
+struct sink {
+  source<T>& _upstream;
+  std::uint64_t _version;
+  std::function<void(T&)> _cb;
+
+  void maybe_update() {
+    if (_version != _upstream._version) {
+      _cb(_upstream.get());
+      _version = _upstream._version;
+    }
+  }
+};
+
+}  // namespace autoopt

lib/importer.cpp

@@ -0,0 +1,73 @@
+#include "autoopt/importer.hpp"
+
+namespace autoopt {
+std::generator<std::pair<double, double>> parse_dat_file(
+    const std::filesystem::path& filename) {
+  // check if file exists
+  if (!std::filesystem::exists(filename)) {
+    throw std::runtime_error("File does not exist: " + filename.string());
+  }
+
+  std::fstream file(filename, std::ios::in);
+
+  std::string line;
+  while (std::getline(file, line)) {
+    if (line.empty() || line == "\r" || line[0] == '#') {
+      continue;
+    }
+    std::istringstream iss(line);
+    double x, y;
+    iss >> x;
+    for (size_t i = 1; i < 5; ++i) {
+      iss >> y;  // skip unused columns
+    }
+
+    co_yield {x, -y};
+  }
+  co_return;
+}
+
+std::pair<std::vector<double>, std::vector<double>> parse_row(auto& iter,
+                                                              auto sentinel) {
+  std::vector<double> xs;
+  std::vector<double> ys;
+  xs.push_back((*iter).first);
+  ys.push_back((*iter).second);
+  ++iter;
+  double last_x = xs[0];
+  bool forward = last_x < (*iter).first;
+  while (iter != sentinel && ((forward && last_x < (*iter).first) ||
+                              (!forward && last_x > (*iter).first))) {
+    xs.push_back((*iter).first);
+    ys.push_back((*iter).second);
+    last_x = (*iter).first;
+    ++iter;
+  }
+  if (!forward) {
+    std::reverse(xs.begin(), xs.end());
+    std::reverse(ys.begin(), ys.end());
+  }
+  return std::make_pair(std::move(xs), std::move(ys));
+}
+
+raw_data read_data(const std::filesystem::path &filename){
+  auto lines = parse_dat_file(filename);
+
+  raw_data result;
+
+  auto it = lines.begin();
+  while (it != lines.end()) {
+    auto [xs, ys] = parse_row(it, lines.end());
+    if (result._xs.rows() == 0) {
+      result._xs = Eigen::VectorXd(xs.size());
+      std::copy(xs.begin(), xs.end(), result._xs.begin());
+    }
+    Eigen::VectorXd ys_vec = Eigen::VectorXd(ys.size());
+    std::copy(ys.begin(), ys.end(), ys_vec.begin());
+    result._yss.push_back(ys_vec);
+  }
+
+  return result;
+}
+
+}  // namespace autoopt

lib/optimizer.cpp

@@ -0,0 +1,224 @@
+#include "autoopt/optimizer.hpp"
+
+#include <Eigen/src/Core/Matrix.h>
+#include <fftw3.h>
+
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+#include <print>
+#include <stdexcept>
+#include <tuple>
+#include <variant>
+
+#include "autoopt/ellipse_new.hpp"
+#include "autoopt/parabola_new.hpp"
+
+namespace autoopt {
+
+void optimizer::import_file(std::filesystem::path file_path) {
+  _data = data{};
+  _data->_raw = read_data(file_path);
+  _data->lb = *_data->_raw._xs.begin();
+  _data->ub = *(_data->_raw._xs.end() - 1);
+  plot_input_cb();
+  init_bounds_cb();
+  update_test_data();
+  update_table_cb();
+}
+
+void optimizer::set_lower_bound(double lb) {
+  if (_data) {
+    _data->lb = lb;
+  }
+  update_test_data();
+}
+
+void optimizer::set_upper_bound(double ub) {
+  if (_data) {
+    _data->ub = ub;
+  }
+  update_test_data();
+}
+
+void optimizer::set_curve_type(int ct) {
+  if (ct == _curve.index()) return;
+  switch (ct) {
+    case CurveType::Parabola:
+      _curve = parabola_new{};
+      break;
+    case CurveType::Ellipse:
+      _curve = ellipse_new{};
+      break;
+  }
+  update_table_cb();
+}
+
+void optimizer::set_initial(size_t i, double value) {
+  _fit_result.reset();
+  std::visit(
+      [&](auto& curve) -> void {
+        if (i >= curve.initial.size()) {
+          throw std::runtime_error("index out of bounds");
+        }
+        curve.initial(i) = value;
+      },
+      _curve);
+  plot_slope_cb();
+}
+
+void optimizer::set_delta(size_t i, double value) {
+  _fit_result.reset();
+  std::visit(
+      [&](auto& curve) -> void {
+        if (i >= curve.delta.size()) {
+          throw std::runtime_error("index out of bounds");
+        }
+        curve.delta(i) = value;
+      },
+      _curve);
+  plot_slope_cb();
+}
+
+void optimizer::flip_curve() {
+  std::visit([](auto& curve) { curve.flip = !curve.flip; }, _curve);
+  plot_slope_cb();
+}
+
+void optimizer::calc_fit() {
+  _fit_result =
+      std::visit([&](auto& curve) { return curve.fit(*_test_data); }, _curve);
+  plot_slope_cb();
+  for (auto param : _fit_result->first) {
+    std::println("param: {}", param);
+  }
+  plot_error_cb();
+}
+
+void optimizer::set_plot_input_cb(
+    std::function<void(Eigen::VectorXd xs, std::vector<Eigen::VectorXd> yss)>
+        cb) {
+  _plot_input_cb = cb;
+}
+
+void optimizer::set_init_bounds_cb(
+    std::function<void(double lb, double ub, double step)> cb) {
+  _init_bounds_cb = cb;
+}
+
+void optimizer::set_plot_slope_cb(
+    std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys,
+                       Eigen::VectorXd ys_curve)>
+        cb) {
+  _plot_slope_cb = cb;
+}
+
+void optimizer::set_update_table_cb(
+    std::function<
+        void(std::vector<std::tuple<std::string, double, double>> entries)>
+        cb) {
+  _param_table_cb = cb;
+}
+
+void optimizer::set_plot_error_cb(
+    std::function<void(Eigen::VectorXd xs, Eigen::VectorXd ys)> cb) {
+  _plot_error_cb = cb;
+}
+
+void optimizer::update_test_data() {
+  _fit_result.reset();
+  if (!_data) {
+    _test_data.reset();
+    return;
+  }
+  auto lb_it = std::lower_bound(_data->_raw._xs.begin(), _data->_raw._xs.end(),
+                                _data->lb);
+  auto ub_it = std::upper_bound(_data->_raw._xs.begin(), _data->_raw._xs.end(),
+                                _data->ub);
+  ssize_t n_samples = std::distance(lb_it, ub_it);
+  if (n_samples < 1) {
+    _test_data.reset();
+    plot_input_cb();
+    return;
+  }
+
+  size_t start_index = std::distance(_data->_raw._xs.begin(), lb_it);
+
+  Eigen::VectorXd xs(n_samples);
+  Eigen::VectorXd ys = Eigen::VectorXd::Zero(n_samples);
+
+  std::copy(lb_it, ub_it, xs.begin());
+  for (size_t i = 0; i < n_samples; i++) {
+    for (auto r : _data->_raw._yss) {
+      ys(i) += r(start_index + i);
+    }
+    ys(i) /= _data->_raw._yss.size();
+  }
+  double center = (*lb_it + *std::prev(ub_it)) / 2;
+  _test_data.emplace(xs, ys, center);
+  plot_slope_cb();
+}
+
+void optimizer::init_bounds_cb() {
+  if (_init_bounds_cb) {
+    _init_bounds_cb(_data->_raw._xs(0), _data->_raw._xs.tail(1)(0),
+                    _data->_raw._xs(1) - _data->_raw._xs(0));
+  }
+}
+
+void optimizer::plot_input_cb() {
+  if (_plot_input_cb) {
+    _plot_input_cb(_data->_raw._xs, _data->_raw._yss);
+  }
+}
+
+void optimizer::plot_slope_cb() {
+  if (_plot_slope_cb) {
+    auto ys_curve =
+        std::visit([&](auto curve) { return curve.plot(*_test_data); }, _curve);
+    _plot_slope_cb(_test_data->_xs, _test_data->_ys, ys_curve);
+  }
+}
+
+void optimizer::update_table_cb() {
+  std::vector<std::tuple<std::string, double, double>> res;
+  std::visit(
+      [&]<typename T>(T& curve) {
+        for (size_t i = 0; i < curve.param_names.size(); i++) {
+          res.emplace_back(curve.param_names[i], curve.initial[i],
+                           curve.delta[i]);
+        }
+      },
+      _curve);
+  _param_table_cb(res);
+}
+
+void optimizer::plot_error_cb() {
+  if (_plot_error_cb) {
+    Eigen::VectorXd ys_curve = std::visit(
+        [&](auto curve) { return curve.plot(*_test_data, _fit_result->first); },
+        _curve);
+    Eigen::VectorXd error = _test_data->_ys - ys_curve;
+    _plot_error_cb(_test_data->_xs, error);
+
+    
+
+    fftw_complex* out = fftw_alloc_complex(error.size());
+    double* input = &error(0);
+    fftw_plan plan = fftw_plan_dft_r2c_1d(error.size(), input, out, FFTW_FORWARD);
+    fftw_execute(plan);
+    fftw_destroy_plan(plan);
+    Eigen::VectorXd out_abs(error.size());
+    for (size_t i = 0; i < error.size(); i++) {
+      out_abs(i) = out[i][0] * out[i][0] + out[i][1]*out[i][1];
+    }
+    fftw_free(out);
+    std::println("begin fftw");
+    for (auto v : out_abs) {
+      std::println("  {}", v);
+    }
+    std::println("end fftw");
+  }
+}
+
+}  // namespace autoopt

pyproject.toml

@@ -0,0 +1,36 @@
+[build-system]
+requires = ["scikit-build-core>=0.11", "nanobind"]
+build-backend = "scikit_build_core.build"
+
+
+[tool.scikit-build]
+# Protect the configuration against future changes in scikit-build-core
+minimum-version = "0.4"
+
+# Setuptools-style build caching in a local directory
+build-dir = "build/{wheel_tag}"
+
+[project]
+name = "slopefit"
+version = "0.1.0"
+authors = [
+  { name="Valentin Leon Stöcker", email="valentin@stoeckermail.de" },
+]
+description = "fit conic sections"
+readme = "README.md"
+requires-python = ">=3.9"
+dependencies = [
+  "numpy",
+  "matplotlib"
+  "pyside6"
+]
+classifiers = [
+    "Programming Language :: Python :: 3",
+    "Operating System :: OS Independent",
+]
+license = "MIT"
+
+[project.urls]
+Homepage = "https://gitea.valentinstoecker.de/vls/slopfit"
+Issues = "https://gitea.valentinstoecker.de/vls/slopfit/issues"
+

src/module.cpp

@@ -0,0 +1,53 @@
+#include <nanobind/eigen/dense.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/stl/filesystem.h>
+#include <nanobind/stl/optional.h>
+#include <nanobind/stl/pair.h>
+#include <nanobind/stl/tuple.h>
+#include <nanobind/stl/string.h>
+#include <nanobind/stl/vector.h>
+#include <nanobind/stl/function.h>
+
+#include <Eigen/Eigen>
+#include <autoopt/optimizer.hpp>
+#include <autoopt/util.hpp>
+#include <format>
+#include <optional>
+
+namespace nb = nanobind;
+
+namespace std {
+
+template <typename T>
+struct formatter<optional<T>> {
+  constexpr auto parse(format_parse_context& ctx) { return ctx.begin(); }
+
+  auto format(const optional<T>& opt, format_context& ctx) const {
+    if (opt) {
+      return format_to(ctx.out(), "{}", opt.value());
+    } else {
+      return format_to(ctx.out(), "None");
+    }
+  }
+};
+
+}  // namespace std
+
+NB_MODULE(slopefit, m) {
+  nb::class_<autoopt::optimizer>(m, "Optimizer")
+      .def(nb::init<>())
+      .def("import_file", &autoopt::optimizer::import_file)
+      .def("set_lower_bound", &autoopt::optimizer::set_lower_bound)
+      .def("set_upper_bound", &autoopt::optimizer::set_upper_bound)
+      .def("set_curve_type", &autoopt::optimizer::set_curve_type)
+      .def("set_initial", &autoopt::optimizer::set_initial)
+      .def("set_delta", &autoopt::optimizer::set_delta)
+      .def("flip_curve", &autoopt::optimizer::flip_curve)
+      .def("calc_fit", &autoopt::optimizer::calc_fit)
+      .def("set_plot_input_cb", &autoopt::optimizer::set_plot_input_cb)
+      .def("set_init_bounds_cb", &autoopt::optimizer::set_init_bounds_cb)
+      .def("set_plot_slope_cb", &autoopt::optimizer::set_plot_slope_cb)
+      .def("set_update_table_cb", &autoopt::optimizer::set_update_table_cb)
+      .def("set_plot_error_cb", &autoopt::optimizer::set_plot_error_cb);
+}

tests/ellipse.cpp

@@ -1,72 +0,0 @@
-#include <gtest/gtest.h>
-
-#include <autoopt/ellipse.hpp>
-#include <autoopt/optimization_problem.hpp>
-#include <autoopt/util.hpp>
-#include <autoopt/btls.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}};
-
-  Eigen::VectorX<double> params(4);
-  params << 100, 1000, deg2rad(1.0), 0.0;
-
-  auto loss_func = [&data_points]<typename T>(const Eigen::VectorX<T>& 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, decltype(loss_func)> problem(loss_func, params);
-
-  auto grad = problem.gradient(params);
-
-  std::cout << grad << std::endl;
-
-  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)
-
-  auto hess = problem.hessian(params);
-
-  std::cout << hess << std::endl;
-
-  Eigen::VectorX<double> params_delta(4);
-  params_delta << 1.0, 1.0, deg2rad(0.1), deg2rad(0.1);
-
-  params(0) += 0.9;  // left_arm
-
-  // log barrier
-  log_barrier_optimization_problem<double> log_barrier_problem(
-      problem,
-      params_delta,
-      1e-3);
-  
-  auto log_barrier_grad = log_barrier_problem.gradient(params);
-
-  std::cout << "Log barrier gradient:" << std::endl;
-  std::cout << log_barrier_grad << std::endl;
-
-  btls(problem);
-}

tests/quadric.cpp

@@ -1,7 +1,6 @@
 #include <gtest/gtest.h>
 
 #include "autoopt/quadric.hpp"
-#include <iostream>
 
 using namespace autoopt;
 

tests/reacitve.cpp

@@ -0,0 +1,21 @@
+#include <gtest/gtest.h>
+
+#include <autoopt/reactive.hpp>
+
+TEST(Reactive, Simple) {
+  using namespace autoopt;
+
+  source<int> s1(0);
+  source<int> s2(0);
+  source<int> s3(0);
+  transform<int, int, int> t1([](int x, int y) { return x + y; }, s1,s2);
+  transform<int, int, int> t2([](int x, int y) { return x + y; }, s2,s3);
+
+  ASSERT_EQ(t2.get(), 0);
+  ASSERT_TRUE(s1.set(1));
+  ASSERT_EQ(t2.get(), 1);
+  s2.set(2);
+  ASSERT_EQ(t2.get(), 3);
+  s3.set(3);
+  ASSERT_EQ(t2.get(), 6);
+}

tests/test_rework.cpp

@@ -0,0 +1,49 @@
+#include <gtest/gtest.h>
+
+#include <autoopt/parabola_new.hpp>
+#include <print>
+#include "autoopt/util.hpp"
+
+
+using namespace autoopt;
+
+
+TEST(New, Parabola) {
+  parabola_new parab;
+  parab.flip = false;
+  parab.initial(0) = 50.0;
+  parab.initial(1) = autoopt::deg2rad(1.0);
+  parab.initial(2) = autoopt::deg2rad(0.05);
+
+  parab.delta(0) = 10.0;
+  parab.delta(1) = autoopt::deg2rad(0.1);
+  parab.delta(2) = autoopt::deg2rad(0.1);
+
+  test_data td;
+  td._xs = Eigen::VectorXd(21);
+  td._ys = Eigen::VectorXd(21);
+  td._center = 10.0;
+
+  for (size_t i = 0; i < 21; i++) {
+    td._xs(i) = static_cast<double>(i);
+  }
+
+  auto res = parab.plot(td);
+
+  for (auto y : res) {
+    std::println("y: {}", y);
+  }
+
+  td._ys = res;
+
+  parab.initial(0) = 54.0;
+  parab.initial(1) = autoopt::deg2rad(1.05);
+
+  auto [fit, rms] = parab.fit(td);
+
+  for (auto y : fit) {
+    std::println("fit: {}", y);
+  }
+
+  std::println("RMS: {}", rms);
+}

todo.md

@@ -0,0 +1,9 @@
+# main class
+
+methods:
+- import_file
+- adjust_bounds
+- select_curve
+- set_param
+- flip
+- calc_fit

ui/__init__.py

@@ -0,0 +1 @@
+

ui/app.py

@@ -0,0 +1,221 @@
+import PySide6.QtWidgets as qw
+import PySide6.QtCore as qc
+from PySide6.QtCore import Slot, Qt
+import PySide6.QtGui as qg
+import pyqtgraph as pg
+
+import sys
+
+sys.path.append('../build')
+
+import slopefit as sf
+
+class ParamTableModel(qc.QAbstractTableModel):
+  HEADERS = ['property', 'initial value', 'delta']
+  init_changed = qc.Signal(int, float)
+  delta_changed = qc.Signal(int, float)
+
+  def __init__(self, parent = None):
+    super().__init__(parent)
+    self._rows = []
+  
+  def rowCount(self, parent=qc.QModelIndex()):
+    return len(self._rows)
+  
+  def columnCount(self, parent=qc.QModelIndex()):
+    return len(self.HEADERS)
+  
+  def data(self, index, role=Qt.ItemDataRole.DisplayRole):
+    if not index.isValid():
+      return None
+    row, col = index.row(), index.column()
+    value = self._rows[row][col]
+
+    if role in (Qt.ItemDataRole.DisplayRole, Qt.ItemDataRole.EditRole):
+      return "" if value is None else str(value)
+    if role == Qt.ItemDataRole.TextAlignmentRole and col > 0:
+      return int(Qt.AlignmentFlag.AlignRight | Qt.AlignmentFlag.AlignVCenter)
+    return None
+  
+  def setData(self, index, value, role=Qt.ItemDataRole.EditRole):
+    if not index.isValid() or role != Qt.ItemDataRole.EditRole:
+      return False
+    row, col = index.row(), index.column()
+    if col == 0:                       # name column: plain text
+      return False
+    else:                              # x / dx columns: float
+      try:
+        self._rows[row][col] = float(value) if value != "" else None
+      except ValueError:
+        return False
+    self.dataChanged.emit(index, index, [Qt.ItemDataRole.DisplayRole, Qt.ItemDataRole.EditRole])
+    if col == 1:
+      self.init_changed.emit(row, float(value))
+    elif col == 2:
+      self.delta_changed.emit(row, float(value))
+    return True
+  
+  def flags(self, index):
+    if not index.isValid() or index.column() == 0:
+      return Qt.ItemFlag.NoItemFlags
+    return Qt.ItemFlag.ItemIsEnabled | Qt.ItemFlag.ItemIsSelectable | Qt.ItemFlag.ItemIsEditable
+
+  def headerData(self, section, orientation, role=Qt.ItemDataRole.DisplayRole):
+    if role == Qt.ItemDataRole.DisplayRole and orientation == Qt.Orientation.Horizontal:
+      return self.HEADERS[section]
+    return None
+
+  def set_params(self, params):
+    self.beginResetModel()
+    self._rows = [[name, v, d] for (name, v, d) in params]
+    self.endResetModel()
+
+
+
+class MainWindow(qw.QMainWindow):
+  def __init__(self):
+    super().__init__()
+    self.setWindowTitle('slopefit')
+
+    
+    self.input_plot = pg.PlotWidget()
+
+    self.lower_bounds = qw.QDoubleSpinBox()
+    self.lower_bounds.setValue(0.0)
+    self.lower_bounds.valueChanged.connect(self.set_lb)
+
+    self.upper_bounds = qw.QDoubleSpinBox()
+    self.upper_bounds.setValue(0.0)
+    self.upper_bounds.valueChanged.connect(self.set_ub)
+
+
+    bounds_layout = qw.QHBoxLayout()
+    bounds_layout.addWidget(self.lower_bounds)
+    bounds_layout.addWidget(self.upper_bounds)
+
+    bounds_widget = qw.QWidget()
+    bounds_widget.setLayout(bounds_layout)
+
+    self.slope_plot = pg.PlotWidget()
+
+    self.curve_box = qw.QComboBox()
+    self.curve_box.addItems(['parabola', 'ellipse'])
+    self.curve_box.currentIndexChanged.connect(self.set_curve_i)
+
+    self.table = qw.QTableView()
+    self.param_model = ParamTableModel()
+    self.param_model.init_changed.connect(self.set_init)
+    self.param_model.delta_changed.connect(self.set_delta)
+    self.table.setModel(self.param_model)
+    self.table.horizontalHeader().setSectionResizeMode(qw.QHeaderView.Stretch)
+
+    self.flip_curve = qw.QCheckBox()
+    self.flip_curve.checkStateChanged.connect(self.set_check)
+
+    self.fit_button = qw.QPushButton()
+    self.fit_button.setText("calculate fit")
+    self.fit_button.clicked.connect(self.do_fit)
+
+    self.error_plot = pg.PlotWidget()
+    plot_group_1 = qw.QHBoxLayout()
+    plot_group_1.addWidget(self.input_plot)
+    plot_group_1.addWidget(self.slope_plot)
+
+    layout = qw.QVBoxLayout()
+    layout.addWidget(qw.QWidget(layout=plot_group_1))
+    layout.addWidget(bounds_widget)
+    layout.addWidget(self.flip_curve)
+    layout.addWidget(self.curve_box)
+    layout.addWidget(self.table)
+    layout.addWidget(self.fit_button)
+    layout.addWidget(self.error_plot)
+
+    central_widget = qw.QWidget()
+    central_widget.setLayout(layout)
+    self.setCentralWidget(central_widget)
+
+    self.init_menu()
+    
+    self.optimizer = sf.Optimizer()
+    self.optimizer.set_plot_input_cb(self.plot_input_cb)
+    self.optimizer.set_init_bounds_cb(self.init_bounds_cb)
+    self.optimizer.set_plot_slope_cb(self.plot_slope_cb)
+    self.optimizer.set_update_table_cb(self.update_table_cb)
+    self.optimizer.set_plot_error_cb(self.plot_error_cb)
+
+  @Slot(float)
+  def set_lb(self, lb):
+    self.optimizer.set_lower_bound(lb)
+    
+  @Slot(float)
+  def set_ub(self, ub):
+    self.optimizer.set_upper_bound(ub)
+  
+  @Slot(int)
+  def set_curve_i(self, i):
+    self.optimizer.set_curve_type(i)
+
+  @Slot(int, float)
+  def set_init(self, i, v):
+    self.optimizer.set_initial(i, v)
+  
+  @Slot(int, float)
+  def set_delta(self, i, v):
+    self.optimizer.set_delta(i, v)
+
+  @Slot(bool)
+  def set_check(self, flip):
+    self.optimizer.flip_curve()
+
+  @Slot()
+  def do_fit(self):
+    self.optimizer.calc_fit()
+
+  def init_menu(self):
+    import_action = qg.QAction('&import', self)
+    import_action.triggered.connect(self.import_file)
+
+    menu = self.menuBar()
+    file_menu = menu.addMenu('&file')
+    file_menu.addAction(import_action)
+
+  def import_file(self):
+    print("import file...")
+    fn, _ = qw.QFileDialog.getOpenFileName(self, caption='import file...', filter='data files (*.dat)')
+    print(f"open file {fn}")
+    self.optimizer.import_file(fn)
+
+  def plot_input_cb(self, xs, yss):
+    self.input_plot.clear()
+    for ys in yss:
+      self.input_plot.plot(xs, ys)
+
+  def init_bounds_cb(self, lb, ub, step):
+    self.lower_bounds.setRange(lb, ub)
+    self.lower_bounds.setSingleStep(step)
+    self.lower_bounds.setValue(lb)
+
+    self.upper_bounds.setRange(lb, ub)
+    self.upper_bounds.setSingleStep(step)
+    self.upper_bounds.setValue(ub)
+  
+  def plot_slope_cb(self, xs, ys_data, ys_curve):
+    self.slope_plot.clear()
+    self.slope_plot.plot(xs, ys_data)
+    self.slope_plot.plot(xs, ys_curve)
+  
+  def update_table_cb(self, entries):
+    self.param_model.set_params(entries)
+    print(entries)
+
+  def plot_error_cb(self, xs, ys):
+    self.error_plot.clear()
+    self.error_plot.plot(xs, ys)
+
+
+app = qw.QApplication(sys.argv)
+
+window = MainWindow()
+window.show()
+
+app.exec()

ui/model.py

@@ -0,0 +1,4 @@
+import PySide6.QtWidgets
+
+
+